C***********************************************************************
C***********************************************************************
      SUBROUTINE WFINRN(NRPLWV,NBANDS,NKPTS,DATAKE,NBANOC,CPTWFP,
     &     NPLWKP,PLWVEN,NFIND,CELEN,NGX, NGY, NGZ, 
     &     LPCTX, LPCTY, LPCTZ, NINDPW, RECC, CPTDUM,CCTDUM)
C=======================================================================
C
C """"""""""""""""""""""""" SUBROUTINE WFINIT """"""""""""""""""""""""""
C
C THIS SUBROUTINES INITIALISES THE WAVEFUNCTIONS FOR THE MOLECULAR
C DYNAMICS BY FILLING THE BANDS AT EACH K POINT WITH THE LOWEST ENERGY
C PLANE WAVE BASIS STATES. IF THE SYSTEM HAS HIGH SYMMETRY THIS METHOD
C MAY NOT BE APPROPRIATE BECAUSE THE LOWEST ENERGY PLANE WAVES MAY NOT
C SPAN THE LOWEST BANDS.
C
C=======================================================================
      IMPLICIT COMPLEX (C)
      EXTERNAL RANF
      DIMENSION DATAKE(*)
      DIMENSION CPTWFP(*),CPTDUM(*),CCTDUM(*)
      DIMENSION CELEN(*)
      DIMENSION LPCTX(*), LPCTY(*), LPCTZ(*), NINDPW(*), RECC(3,3)
C=======================================================================
C
C                       DIMENSION STATEMENTS
C
C PLWVEN = THE ENERGIES OF THE LOWEST KINETIC ENERGY PLANE WAVE BASIS
C          STATES
C NFIND = THE INDICES OF THE LOWEST ENERGY PLANE WAVE BASIS STATES
C=======================================================================
      DIMENSION PLWVEN(*)
      DIMENSION NFIND(*)
      DATA HSQDTM /3.810033/
      DATA ISEED /-1/
C=======================================================================
C INITIALISE THE ENERGIES OF THE BASIS STATES TO A LARGE VALUE SO THAT
C LOWER ENERGY BASIS STATES WILL BE FOUND
C=======================================================================
      WRITE (*,*) 'WELCOME TO WFINIT-RANDOM'
      NBASE = NBANDS * 3
      DO 2010 MM = 1,NBASE
        PLWVEN(MM) = 10000.0
 2010 CONTINUE
C=======================================================================
C RUN THROUGH THE SET OF PLANE WAVE BASIS STATES AT THE PRESENT K POINT
C AND FIND THE ENERGY .
C=======================================================================
      DO 2100 NNN = 1,NPLWKP
        NINDX = 1 + ( 7 * (NNN - 1) )
C======================================================================
C IF THE ENERGY OF THIS BASIS STATES IS HIGHER THAN THE BANDS LOWEST 
C ENERGIES FOUND SO FAR MOVE ONTO THE NEXT BASIS STATES
C======================================================================
        IF (DATAKE(NINDX).GT.PLWVEN(NBASE)) GOTO 2108
C======================================================================
C OTHERWISE RUN THROUGH ALL PLANE WAVE BASIS FROM WHICH THE INITIALIZED
C WAVEFUNCTION IS BEING GENERATED
C=====================================================================
        DO 2110 M = 1,NBASE
          IF (PLWVEN(M).GT.DATAKE(NINDX)) GO TO 2200
 2110   CONTINUE
        GO TO 2108
 2200   NSWAP = M
C=======================================================================
C PUT THIS BASIS STATE AT ITS RIGHTFULL PLACE AMONGST THE LOWEST ENERGY
C PLANE WAVE BASIS STATES AND MOVE THE HIGHER ENERGY COMPONENTS UP
C ONE POSITION
C======================================================================
        DO 2109 M = NBASE,NSWAP+1,-1
          IF (M.GT.NBASE) GO TO 2109
          PLWVEN(M) = PLWVEN(M-1)
          NFIND(M) = NFIND(M-1)
 2109   CONTINUE
        PLWVEN(NSWAP) = DATAKE(NINDX)
        NFIND(NSWAP) = NNN
 2108   CONTINUE
 2100 CONTINUE
C======================================================================
C     The plane waves which are to be included in the initial
C     wave functions. Find their GX, GY, GZ 
C======================================================================
C
      DO 2450 NN = 1, NBASE
        M = NINDPW(NFIND(NN))
        IGX = LPCTX((M-1)       - ( (M-1) /   NGX     ) * NGX + 1)
        IGY = LPCTY((M-1) / NGX - ( (M-1) / (NGX*NGY) ) * NGY + 1)
        IGZ = LPCTZ((M-1) / (NGX * NGY)                       + 1)
C        WRITE (*,2502) IGX, IGY, IGZ, PLWVEN(NN)
 2450 CONTINUE
 2502 FORMAT(1X,'  INITIAL WAVEFUNCTION INCLUDES: G =', 3I4, F15.6)
C=======================================================================
C FROM THE INDICES OF THE LOWEST ENERGY BASIS STATES INITIALISE THE
C WAVEFUNCTIONS
C=======================================================================
C
C      RANDOM INTIALISATION IS ADDED HERE
C
C-----------------------------------------------------------------------
C
C
      DO 2500 NN = 1,NBANDS
        NINDX = NRPLWV*(NN-1)
        DO 2500 NN1 = 1,NBASE
          IF (PLWVEN(NN1).EQ.0.0) THEN
            PLWVEN(NN1) = 1.0
          ENDIF
          CPTWFP(NFIND(NN1)+NINDX) = (0.5,0.0) *
     &                    RANF(ISEED) / (PLWVEN(NN1)**2)
2500  CONTINUE
C======================================================================
C
C     NORMALIZE THE FIRST BAND IN ORDER TO ORTHOGONAL THE HIGHER BANDS
C
C======================================================================
      DO 2600 NN = 1,NBANDS
        NINDW = NRPLWV * ( NN - 1 )
        DO 2601 M = 1,NPLWKP
          CPTDUM(M) = CPTWFP(M+NINDW)
2601    CONTINUE
        DO 2602 M = 1,NPLWKP
          CCTDUM(M) = CONJG(CPTDUM(M))
2602    CONTINUE
        CWFMAG = (0.0,0.0)
        DO 2603 M = 1,NPLWKP
          CWFMAG = CWFMAG + CCTDUM(M) * CPTDUM(M)
2603    CONTINUE
        WFMAG = REAL(CWFMAG)
        IF (WFMAG.GT.0) GO TO 2605
        WRITE (6,2604) NN,M
        STOP
2605  CONTINUE
        WFMINV = 1.0 / SQRT(WFMAG)
        DO 2606 M = 1,NPLWKP
          CPTDUM(M) = WFMINV * CPTDUM(M)
          CCTDUM(M) = CONJG(CPTDUM(M))
2606    CONTINUE
        DO 2607 M = 1,NPLWKP
          CPTWFP(M+NINDW) = CPTDUM(M)
2607    CONTINUE
C======================================================================
C
C     NOW ORTHOGONALIZE THE WAVEFUNCTION OF ALL THE HIGHER BAND TO THE
C     PRESENT NORMALIZED WAVEFUNCTION
C
C======================================================================
        DO 2608 MM = NN+1,NBANDS
          NINDX = NRPLWV * ( MM - 1 )
          CSCPD = (0.0,0.0)
          DO 2609 M = 1,NPLWKP
            CSCPD = CSCPD + CCTDUM(M) * CPTWFP(M+NINDX)
2609      CONTINUE
          DO 2610 M = 1,NPLWKP
            CPTWFP(M+NINDX) = CPTWFP(M+NINDX) - CSCPD * CPTDUM(M)
2610      CONTINUE
2608    CONTINUE
2600  CONTINUE
2604  FORMAT(1X,'WAVEFUNCTION LINEARLY DEPENDENT',2I5)
C======================================================================
C
C      We now calculate the kinetic energy of each band
C
C======================================================================
      DO 2800 NN = 1,NBANDS 
        NINDX = NRPLWV * ( NN - 1 )
        CELEN(NN) = (0.0,0.0)
        DO 2810 M = 1,NPLWKP
          CELEN(NN) = CELEN(NN) +
     &    DATAKE(1+7*(M-1)) * CPTWFP(M+NINDX) * CONJG(CPTWFP(M+NINDX))
 2810   CONTINUE
C        WRITE(*,*)'  KINETIC ENERGY IS: ', CELEN(NN)
 2800 CONTINUE
      RETURN
      END


      FUNCTION RANF(IDUM)
      PARAMETER (M=714025,IA=1366,IC=150889,RM=1.0/M)
      DIMENSION IR(97)
      DATA IFF /0/
      IF(IDUM.LT.0.OR.IFF.EQ.0) THEN
        IFF=1
        IDUM=MOD(IC-IDUM,M)
        DO 11 J=1,97
           IDUM=MOD(IA*IDUM+IC,M)
           IR(J)=IDUM
  11    CONTINUE
        IDUM=MOD(IA*IDUM+IC,M)
        IY=IDUM
      ENDIF
      J=1+(97*IY)/M
      IF(J.GT.97.OR.J.LT.1) THEN
      PRINT *,' RANF: J=',J
      STOP
      ENDIF
      IY=IR(J)
      RANF=IY*RM
      IDUM=MOD(IA*IDUM+IC,M)
      IR(J)=IDUM
      RETURN
      END
