      PROGRAM GRID
C
C Create kpoint grid in CASTEP or CETEP style.
C Optionally remove symmetry related points
C
C Output file: kpts.out
C
C Compile with -r8
C
C MJR 20/3/97
C
      IMPLICIT NONE
      INTEGER NMAX
      PARAMETER (NMAX=1000)
      INTEGER NGX,NGY,NGZ,NGZMAX,I,J,K
      REAL DX,DY,DZ,KPT(NMAX,4)
      INTEGER NX,NY,NZ,INV,NKPTS
      REAL X,Y,Z,WT,DUMMY
      LOGICAL CETEP
C
      WRITE(*,*)' WELCOME TO K-POINT GENERATOR'
      WRITE(*,*)' PLEASE INPUT REQUIRED GRID SIZE'
      WRITE(*,*)' NGX, NGY, NGZ'
      READ(*,*) NGX,NGY,NGZ
      NKPTS=NGX*NGY*NGZ
      IF (NKPTS.GT.NMAX) STOP 'INCREASE NMAX'
      WRITE(*,*)
      WRITE(*,*)' PLEASE INPUT OFFSET: ZERO INCLUDES GAMMA POINT'
      WRITE(*,*)'                      0.5 SHIFTS BY 0.5*KPT SPACING'
      WRITE(*,*)' DX,DY,DZ'
      READ(*,*) DX,DY,DZ
      WRITE(*,*)' ASSUME INVERSION (1=YES, 0=NO)'
      READ(*,*) INV
      WRITE(*,*)' CETEP OR CASTEP STYLE OUTPUT? (1=CETEP, 0=CASTEP)'
      CETEP=.TRUE.
      READ(*,*)DUMMY
      IF (DUMMY.EQ.0) CETEP=.FALSE.
C
C GENERATE ALL KPOINTS
C
      WT=1.0/NKPTS
      I=0
      DO NX=0,NGX-1
        DO NY=0,NGY-1
          DO NZ=0,NGZ-1
            I=I+1
            KPT(I,1)=(NX+DX)/NGX
            KPT(I,2)=(NY+DY)/NGY
            KPT(I,3)=(NZ+DZ)/NGZ
            KPT(I,4)=WT
          END DO
        END DO
      END DO
C
C REMOVE INVERSION RELATED KPOINTS
C
      IF (INV.EQ.1) THEN
        DO I=1,NKPTS
          DO J=I+1,NKPTS
            IF (ABS(KPT(I,1)-MOD(2.0-KPT(J,1),1.0)).LT.1E-8
     &        .AND.ABS(KPT(I,2)-MOD(2.0-KPT(J,2),1.0)).LT.1E-8
     &        .AND.ABS(KPT(I,3)-MOD(2.0-KPT(J,3),1.0)).LT.1E-8) THEN
              KPT(I,4)=KPT(I,4)+KPT(J,4)
              KPT(J,4)=0
            ENDIF
          END DO
         END DO
       END IF
C
C WRITE OUT
C
       OPEN(10,FILE='kpts.out')
       J=0
       DO I=1,NKPTS
         IF (KPT(I,4).NE.0) THEN
           IF (CETEP) THEN
              WRITE (10,1000) (KPT(I,K),K=1,4)
           ELSE
              WRITE (10,1000) (KPT(I,K),K=1,3)
           END IF
           J=J+1
         END IF
       END DO
1000   FORMAT(1X,4F15.10)
C
       IF (.NOT.CETEP) THEN
         DO I=1,NKPTS
           IF (KPT(I,4).NE.0) WRITE (10,1000) KPT(I,4)
         END DO
       END IF
C
       WRITE(*,*)'NKPTS=',J
C
       END
