C***********************************************************************
C FILE: GENBTR.FOR
C
C HISTORY: 19/2/90: COPIED FROM GENSP1.FOR
C                   RECC,RECI IMPLEMENTED BY BH
C          25/2/90: DATAKE(2,) (3,) AND (4,) NOW ENERGIES
C                   OFF-DIAGONAL ELEMENTS OF STRESS TENSOR
C          11/3/90: PRINT-MESSAGE ABOUT CUT-OFF
C          12/3/90:         -//-        HOW MANY PLANE-WAVES ACCEPTED 
C          25/5/90: REDUNDANT CAR-PARRINELLO PART REMOVED BY BH
C        05-JUN-90: DNLG AND DNLKG ADDED FOR THE USE OF NON-LOCAL
C		    PSEUDOPOTENTIAL BY X WENG
C        09-JUL-90: DNLG AND DNLKG ARE ALL IN CARTESIAN COORDINATION 
C                   BY X WENG
C***********************************************************************
      SUBROUTINE GENBTR(NRPLWV,NGX,NGY,NGZ,NKPTS,ENMAX,NINDPW,NPLWKP,
     &      VKPT,LPCTX,LPCTY,LPCTZ,DATAKE,RECC,RECI,IPRINT,
     &      DNLG,DNLKG)
C=======================================================================
C
C ]]]]]]]]]]]]]]]]]]]]]] SUBROUTINE GENSP [[[[[[[[[[[[[[[[[[[[[[[[[[[[[[
C
C SUBROUTINE GENSP PERFORMS A NUMBER OF TASKS. THE DATA ARRAYS COSDAT &
C DECDAT FOR THE ELECTRON DYNAMICS ARE INITIALISED. THE INDEXING SYSTEM
C FOR PADDING THE SPHERES OF PLANE WAVES AT EACH K POINT INTO THE BOX
C USED FOR THE FAST FOURIER TRANSFORMS IS ALSO COMPUTED AS ARE THE
C KINETIC ENERGIES OF THE PLANE WAVE BASIS STATES AT EACH K POINT AND
C THE X,Y AND Z COMPONENTS OF THE KINETIC ENERGY. THE INDEXING SYSTEM
C IS BASED ON THE INITIAL SIZE OF THE UNIT CELL THOUGH THE KINETIC
C ENERGIES OF THE PLANE WAVE BASIS STATES OBVIOUSLY CHANGES AS THE
C SIZE OF THE CELL CHANGES
C=======================================================================
      IMPLICIT COMPLEX (C)
      DIMENSION NINDPW(*)
      DIMENSION NPLWKP(*)
      DIMENSION VKPT(3,*)
      DIMENSION LPCTX(*),LPCTY(*),LPCTZ(*)
      DIMENSION DATAKE(*)
      DIMENSION RECC(3,*)
      DIMENSION RECI(3,*)
      DIMENSION DNLG(NRPLWV,3,NKPTS),DNLKG(NRPLWV,0:3,NKPTS)
C=======================================================================
C 	TPI = (NOT SURPRISINGLY) 2*PI
C 	HSQDTM = (PLANCK'S CONSTANT/(2*PI))**2/(2*ELECTRON MASS)
C======================================================================
      DATA TPI,HSQDTM /6.283185307,3.810033/
      TESTMX=0.0
C=======================================================================
C THE FINAL SECTION OF THE SUBROUTINE CALCULATES THE INDEX IN THE
C RECIPROCAL LATTICE GRID OF EACH PLANE WAVE BASIS STATES AT EACH
C K POINT. THE KINETIC ENERGIES OF THE BASIS STATES AND THE X,Y AND Z
C COMPONENTS OF THE KINETIC ENERGY ARE CALCULATED. IF THE CUT-OFF ENERGY
C ENMAX IS TOO LARGE SO THAT THERE ARE MORE THAN NRPLWV PLANE WAVE BASIS
C STATES OF ENERGY LESS THAN ENMAX AT ANY K POINT, ENMAX IS
C AUTOMATICALLY REDUCED
C=======================================================================
 10   CONTINUE
      ACCMXC=0.0
      ACCMXI=0.0
      DO 2000 N=1,NKPTS
        NLBOXI=1
        NSBOXI=1
        DO 2010 NZ=1,NGZ
          GZX=RECC(3,1)*(LPCTZ(NZ)+VKPT(3,N))
          GZY=RECC(3,2)*(LPCTZ(NZ)+VKPT(3,N))
          GZZ=RECC(3,3)*(LPCTZ(NZ)+VKPT(3,N))
          GIZX=RECI(3,1)*(LPCTZ(NZ)+VKPT(3,N))
          GIZY=RECI(3,2)*(LPCTZ(NZ)+VKPT(3,N))
          GIZZ=RECI(3,3)*(LPCTZ(NZ)+VKPT(3,N))
          DO 2011 NY=1,NGY
            GYX=RECC(2,1)*(LPCTY(NY)+VKPT(2,N))
            GYY=RECC(2,2)*(LPCTY(NY)+VKPT(2,N))
            GYZ=RECC(2,3)*(LPCTY(NY)+VKPT(2,N))
            GIYX=RECI(2,1)*(LPCTY(NY)+VKPT(2,N))
            GIYY=RECI(2,2)*(LPCTY(NY)+VKPT(2,N))
            GIYZ=RECI(2,3)*(LPCTY(NY)+VKPT(2,N))
            DO 2012 NX=1,NGX
              GXX=RECC(1,1)*(LPCTX(NX)+VKPT(1,N))
              GXY=RECC(1,2)*(LPCTX(NX)+VKPT(1,N))
              GXZ=RECC(1,3)*(LPCTX(NX)+VKPT(1,N))
              GX=GXX+GYX+GZX
              GY=GXY+GYY+GZY
              GZ=GXZ+GYZ+GZZ
              GIXX=RECI(1,1)*(LPCTX(NX)+VKPT(1,N))
              GIXY=RECI(1,2)*(LPCTX(NX)+VKPT(1,N))
              G1=GX
              G2=GY
              G3=GZ
              GIXZ=RECI(1,3)*(LPCTX(NX)+VKPT(1,N))
              ENERGX=HSQDTM*(GX**2)
              ENERGY=HSQDTM*(GY**2)
              ENERGZ=HSQDTM*(GZ**2)
              ENERG=ENERGX+ENERGY+ENERGZ
              ENERIX=HSQDTM*((GIXX+GIYX+GIZX)**2)
              ENERIY=HSQDTM*((GIXY+GIYY+GIZY)**2)
              ENERIZ=HSQDTM*((GIXZ+GIYZ+GIZZ)**2)
              ENERGI=ENERIX+ENERIY+ENERIZ
              IF(ENERGI.GT.TESTMX) TESTMX=ENERGI
C=======================================================================
C CHECK TO SEE IF THE KINETIC ENERGY OF THE PLANE WAVE IS LESS THAN
C ENMAX IN WHICH CASE THE PLANE WAVE IS INCLUDED IN THE SET OF BASIS
C STATES FOR THE K POINT
C=======================================================================
              IF(ENERGI.GE.ENMAX) GOTO 2017
                IF(ENERG.GT.ACCMXC) ACCMXC=ENERG
                IF(ENERGI.GT.ACCMXI) ACCMXI=ENERGI
                NINDPW(NSBOXI+((N-1)*NRPLWV))=NLBOXI
                NINDX=7*(NSBOXI-1+((N-1)*NRPLWV))
                DATAKE(1+NINDX)=ENERG
                DATAKE(2+NINDX)=HSQDTM*(G1**2)
                DATAKE(3+NINDX)=HSQDTM*(G2**2)
                DATAKE(4+NINDX)=HSQDTM*(G3**2)
                DATAKE(5+NINDX)=HSQDTM*G2*G3
                DATAKE(6+NINDX)=HSQDTM*G3*G1
                DATAKE(7+NINDX)=HSQDTM*G1*G2
C=======================================================================
C     FOR THE USE OF NON-LOCAL PSEUDOPOTENTIAL ONLY
C     DNLG(*,3,*) = (Gx,Gy, Gz) in cartesian co-ordination
C=======================================================================
C was            DNLG(1,NSBOXI,N)=LPCTY(NX)
C                DNLG(2,NSBOXI,N)=LPCTY(NY)    << in unit of reciprocal 
C                DNLG(3,NSBOXI,N)=LPCTZ(NZ)       vectors
C
                 DNLG(NSBOXI,1,N)=GX
                 DNLG(NSBOXI,2,N)=GY
                 DNLG(NSBOXI,3,N)=GZ
C=======================================================================
C   THIS IS IN CARTESIAN COODINATION, TOO.
C=======================================================================
                DUM=SQRT(GX**2+GY**2+GZ**2)
                DNLKG(NSBOXI,0,N)=DUM
                IF(DUM.LE.0.000001) THEN
                  DO 2015 I=1,3
                    DNLKG(NSBOXI,I,N)=0.57735027
 2015             CONTINUE
                ELSE 
                  DNLKG(NSBOXI,1,N)=GX/DUM
                  DNLKG(NSBOXI,2,N)=GY/DUM
                  DNLKG(NSBOXI,3,N)=GZ/DUM
                END IF
                NSBOXI=NSBOXI+1
 2017         CONTINUE
              NLBOXI=NLBOXI+1
 2012       CONTINUE
 2011     CONTINUE
 2010   CONTINUE
C=======================================================================
C CHECK TO SEE IF THERE ARE LESS THAN NRPLWV BASIS STATES AT THIS K POIN
C IF NOT REDUCE ENMAX AND START THE WHOLE PROCESS AGAIN
C=======================================================================
        IF(NRPLWV.LT.NSBOXI) THEN
          WRITE(*,2020) ENMAX,TESTMX,NSBOXI,NRPLWV
          ENMAX=ENMAX-1
          GOTO 10
        ENDIF
        NPLWKP(N)=NSBOXI-1
        IF(IPRINT.GE.1) WRITE(*,2050) (NSBOXI-1),N
 2000 CONTINUE
      WRITE(*,2030) ACCMXI
      WRITE(*,2040) ACCMXC
C
 2020 FORMAT(/,1X,'GENTR: ENMAX,TESTMAX,NSBOXI,NRPLWV',2F12.6,2I7)
 2030 FORMAT(2X,'PLANE-WAVES UP TO ',F8.2,' eV IN ORIGINAL CELL ',
     &            'HAVE BEEN ACCEPTED ')
 2040 FORMAT(2X,'PLANE-WAVES UP TO ',F8.2,' eV IN   N E W  CELL ',
     &            'HAVE BEEN ACCEPTED')
 2050 FORMAT(/,' GENTR:',I5,' PLANE-WAVES FOR K-POINT No',I3)
      RETURN
      END
