cryapi_inp.f

      Module Numbers
      Implicit None
      !  A module to look after the data types used by Crystal.
      Integer, Parameter :: float = Selected_real_kind( 13, 100 )
      !  Floating Cardinals 
      End Module Numbers
      
      Module text_module
      INTEGER, PARAMETER, PUBLIC :: NANGECP = 5
      CHARACTER(LEN=24)::  NSPIN(2)
      CHARACTER(LEN=24)::  MSPIN(2)
      CHARACTER(LEN=4) :: NDN(6)
      CHARACTER(LEN=6) :: XMOT(NANGECP) 
      character(len=1) :: tipo(0:1)
      CHARACTER(LEN=24) :: HFTYPE(0:2)
      CHARACTER(LEN=23):: NOMC(8),NOMX(5)
      CHARACTER(LEN=2) :: SYMBAT(0:93)
      DATA NSPIN/'    ALPHA+BETA ELECTRONS','    ALPHA-BETA ELECTRONS'/
      DATA MSPIN/'    ALPHA      ELECTRONS','    BETA       ELECTRONS'/
      DATA NDN/' S',' SP',' P',' D',' F',' G'/
      DATA XMOT/ 'W0 TMS','P0 TMS','P1 TMS','P2 TMS','P3 TMS'/
      DATA TIPO/'C','R'/
      DATA HFTYPE/' RESTRICTED CLOSED SHELL',
     *' RESTRICTED OPEN SHELL  ',' UNRESTRICTED OPEN SHELL'/
      DATA NOMX/'DIRAC-SLATER LDA','VON BARTH-HEDIN','BECKE',
     *'PERDEW-BURKE-ERNZERHOF','PERDEW-WANG GGA'/
      DATA NOMC/'PERDEW-WANG LSD','PERDEW-ZUNGER','VOSKO-WILK-NUSAIR',
     *'VON BARTH-HEDIN','PERDEW 86','PERDEW-WANG GGA','LEE-YANG-PARR',
     *'PERDEW-BURKE-ERNZERHOF'/
      DATA SYMBAT/   'XX','H ','HE','LI','BE','B ','C ','N ',
     1'O ','F ','NE','NA','MG','AL','SI','P ','S ','CL','AR',
     2'K ','CA','SC','TI','V ','CR','MN','FE','CO','NI','CU',
     3'ZN','GA','GE','AS','SE','BR','KR','RB','SR','Y ','ZR',
     4'NB','MO','TC','RU','RH','PD','AG','CD','IN','SN','SB',
     5'TE','I ','XE','CS','BA','LA','CE','PR','ND','PM','SM',
     6'EU','GD','TB','DY','HO','ER','TM','YB','LU','HF','TA',
     7'W ','RE','OS','IR','PT','AU','HG','TL','PB','BI','PO',
     8'AT','RN','FR','RA','AC','TH','PA','U ','QC'/
      end  Module text_module

      Module parame_module
      Use numbers
      Implicit real(float) ( a-h, o-z )
      PARAMETER (LIM001=868,LIM002=15000,LIM015=8000,LIM016=2000,
     *LIM007=3500,LIM006=LIM007+LIM007-1,
     ,LIM018=LIM002*LIM001*2,LIM042=lim002*2)
      End Module parame_module

      Module parinf_module
      Use numbers

      Implicit None

      Integer, Parameter, Public :: limprn = 300
      Integer, Parameter, Public :: limftn = 200

      Real( float ), Dimension( : )       , allocatable,public :: par
      Integer      , Dimension( 1:limprn ), Public :: lprint
      Integer      , Dimension( 1:limftn ), Public :: iunit
      Integer      , Dimension( : )       , allocatable,public :: itol
      Integer      , Dimension( : )       , allocatable,public :: inf
      Character( Len = 80 )               , Public :: itit
      Integer                             , Public :: iin,iout
      End Module parinf_module

      Module basato_module
      Use numbers
      Use parame_module
      Implicit real( float ) ( a-h, o-z )
      COMMON/BASATO/
     * AZNUC(LIM016),XA(3,LIM016),
     * CHE(LIM015),EXAD(LIM015+1),XL(3,LIM015),
     * EXX(LIM042),C1(LIM042),C2(LIM042),C3(LIM042),
     * CMAX(LIM042),C2W(LIM042),C3W(LIM042),
     * NAT(LIM016),NSHPRI(LIM016+1),IPSEUD(LIM016),
     * LAA(LIM015+1),LAN(LIM015),LAT(LIM015),LATAO(LIM015),
     * NDQ(LIM015+1),LATOAT(LIM015)!IASYMMET(LIM016),IPRIM(LIM016)
      End Module basato_module

      MODULE INFPOT_MODULE
        USE NUMBERS
        USE PARAME_MODULE
        INTEGER,PARAMETER :: NANG=5
        REAL(FLOAT),DIMENSION(LIM042) :: APOT,CPOT
        INTEGER,DIMENSION(LIM042) :: NPOT
        INTEGER,DIMENSION(LIM016*NANG) :: NBTYP
        INTEGER,DIMENSION(LIM016+1) :: NSOM
      END MODULE INFPOT_MODULE

      MODULE GVECT_MODULE
        USE NUMBERS
        USE PARAME_MODULE
        REAL(FLOAT),DIMENSION(3,3) :: PARET(3,3),W1R(3,3)
        REAL(FLOAT),DIMENSION(LIM007+1) :: GMODUS
        REAL(FLOAT),DIMENSION(3,LIM006) :: XG
        INTEGER,DIMENSION(LIM007+1) :: NM,MN
        INTEGER,DIMENSION(LIM006) :: NN1
        INTEGER,DIMENSION(3,LIM006) :: LG
      END MODULE GVECT_MODULE

      Module retic_module
      Use numbers
      Use parame_module
      REAL(FLOAT),DIMENSION(3,3), public :: BRET
      REAL(FLOAT),DIMENSION(LIM001), public :: WPJ
      INTEGER,DIMENSION(LIM001), public :: LATVRS
      INTEGER,DIMENSION(3,LIM001), public :: JJ
      INTEGER,DIMENSION(3,3,48), public :: IRR
      INTEGER, public :: IS,ISP,NKF,NKIF,IS1,IS2,IS3
      End Module retic_module

      MODULE EPESI_MODULE
      USE NUMBERS
      USE PARAME_MODULE
      REAL(FLOAT),DIMENSION(LIM018) :: ENE,ALFA
      END MODULE EPESI_MODULE
      
      MODULE MOLBAR_MODULE
      USE NUMBERS
      USE PARAME_MODULE
      REAL(FLOAT),DIMENSION(3,LIM016), public :: XBAR
      INTEGER,DIMENSION(LIM016), public :: N1MOL
      END MODULE MOLBAR_MODULE
      
      Module xyvdim_module
      Use numbers
      Use parame_module
      REAL(FLOAT),DIMENSION(3,3,48), public :: XYV
      REAL(FLOAT),DIMENSION(3,48), public :: TRASV
      INTEGER,DIMENSION(48), public :: NINV
      INTEGER,DIMENSION(48,48), public :: MULTAB
      End Module xyvdim_module
!
!     MODULE ROTMATRIX
!     USE NUMBERS
! ALLOCATABLE VECTORS FOR ROTATION MATRICES (DIMENSION LIMROT1)
!     INTEGER :: LIMROT1
!     REAL(FLOAT),ALLOCATABLE :: TTO(:)
!     INTEGER, ALLOCATABLE :: MMO(:),MMOM(:)
! OTHER RELATED VECTORS & MATRICES (LMAX_DFT7=7)
!     INTEGER :: MINZ(48*(7+1)+1)
!     END MODULE ROTMATRIX

      SUBROUTINE CRYREAD
      USE NUMBERS
      USE PARAME_MODULE
      USE PARINF_MODULE
      USE BASATO_MODULE
      USE INFPOT_MODULE
      USE GVECT_MODULE
      USE XYVDIM_MODULE
      USE TEXT_MODULE
      USE MOLBAR_MODULE
      IMPLICIT REAL(FLOAT) (A-H,O-Z)
      LOGICAL IPRAT(0:399), EXIST
      CHARACTER(LEN=6) :: ZNAMZ='CRYREAD'
      INTEGER,DIMENSION(:),ALLOCATABLE :: IKY
      INTEGER,DIMENSION(:),ALLOCATABLE :: ncf,inzvlb,idime,idimf,idmcou
      INTEGER,DIMENSION(:),ALLOCATABLE :: lav,mgnav,la3,la4,jncdu,irof,
     *jrof
      INTEGER,DIMENSION(:),ALLOCATABLE ::ICCAT,ICCS3,KVRSP,iccv,icct,icc
      INTEGER,DIMENSION(:),ALLOCATABLE :: ina12,nlana
      INTEGER,DIMENSION(:),ALLOCATABLE :: ila12t,jpoint,iccs1,nla21t
      INTEGER,DIMENSION(:),ALLOCATABLE :: nnnc,nnnc2,la34x,la34v,ilana
      INTEGER,DIMENSION(:),ALLOCATABLE :: nshg,ngshg,nqgshg,nstatg
      INTEGER,DIMENSION(:),ALLOCATABLE :: nngi,nshgi
      REAL(FLOAT),DIMENSION(:),ALLOCATABLE :: qtot,sg,fg,pg,a
      IIN=5
      IOUT=6
      iuni=97
c
      INQUIRE(FILE='GRED.DAT',EXIST=EXIST)
      IF(.NOT.EXIST)THEN
         CALL ERRNIC(0,IUNI,ZNAMZ,
     *'FILE GRED.DAT NOT FOUND - NO WF DATA')
         ELSE
           OPEN(IUNI,FILE='GRED.DAT',FORM='FORMATTED',
     *STATUS='OLD')
         ENDIF
      do i=1,limftn
      iunit(i)=i
      enddo
      READ(IUNI,107)ITIT
      READ(IUNI,106)LUMINF,LUMTOL,LUMPAR
      ALLOCATE(INF(1:LUMINF),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of INF')
      ALLOCATE(ITOL(1:LUMTOL),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of ITOL')
      ALLOCATE(PAR(1:LUMPAR),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of PAR')
      CALL IREADF(IUNI,INF,LUMINF)
      CALL IREADF(IUNI,ITOL,LUMTOL)
      CALL RREADF(IUNI,PAR,LUMPAR)
! READ TYPE OF HAMILTONIAN
      WRITE(IOUT,203)ITIT,HFTYPE(INF(27)+INF(64))
      IF(INF(170).NE.0)THEN
      WRITE(IOUT,1007)
      IF((INF(172).NE.0).AND.(INF(171).NE.0))THEN
      WRITE(IOUT,9203)TRIM(NOMX(INF(172))),TRIM(NOMC(INF(171)))
      IF(ABS(PAR(39)-1._FLOAT).GE.1E-11_FLOAT)WRITE(IOUT,2011)PAR(39)
      IF(ABS(PAR(38)-1._FLOAT).GE.1E-11_FLOAT)WRITE(IOUT,2012)PAR(38)
      ELSE
      IF(INF(172).NE.0)THEN
      WRITE(IOUT,9201)TRIM(NOMX(INF(172)))
      IF(ABS(PAR(39)-1._FLOAT).GE.1E-11_FLOAT)WRITE(IOUT,2011)PAR(39)
      ENDIF
      IF(INF(171).NE.0)THEN
      WRITE(IOUT,9202)TRIM(NOMC(INF(171)))
      IF(ABS(PAR(38)-1._FLOAT).GE.1E-11_FLOAT)WRITE(IOUT,2012)PAR(38)
      ENDIF
      ENDIF
      IF(INF(173).NE.0)WRITE(IOUT,1445)PAR(48)
      ELSE
      WRITE(IOUT,1008)
      ENDIF
1007  FORMAT(' KOHN-SHAM HAMILTONIAN')
1008  FORMAT(' HARTREE-FOCK HAMILTONIAN')
9201  FORMAT(/' THE EXCHANGE FUNCTIONAL    ',A,' IS ACTIVE')
9202  FORMAT(/' THE CORRELATION FUNCTIONAL ',A,' IS ACTIVE')
9203  FORMAT(/' (EXCHANGE)[CORRELATION] FUNCTIONAL:(',A,')[',A,']'/)
1445  FORMAT(/' HYBRID EXCHANGE - PERCENTAGE OF FOCK EXCHANGE =',
     *T48,F10.4)
2011  FORMAT(' NON-LOCAL WEIGHTING FACTOR (EXCHANGE) =',T48,F10.4)
2012  FORMAT(' NON-LOCAL WEIGHTING FACTOR [CORRELATION] =',T48,F10.4)
! CHECK OF DIMENSIONS OF STATIC ARRAYS
      LAF=INF(20)
      IF(INF(20).GT.LIM015)CALL ERRNIC(0,INF(20),ZNAMZ,
     *'TOO MANY SHELLS - INCREASE LIM015 TO')
      LAF1=LAF+1
      LAF3=3*LAF
      NAF=INF(24)
      IF(INF(24).GT.LIM016)CALL ERRNIC(0,INF(24),ZNAMZ,
     *'TOO MANY ATOMS - INCREASE LIM016 TO')
      NAF1=NAF+1
      NAF3=3*NAF
      NDF=INF(7)
      IF(INF(7).GT.LIM002)CALL ERRNIC(0,INF(7),ZNAMZ,
     *'TOO MANY AO - INCREASE LIM002 TO')
      NPRIM=INF(75)
      IF(INF(75).GT.LIM042)CALL ERRNIC(0,INF(75),ZNAMZ,
     *'TOO MANY PRIMTIVES - INCREASE LIM042 TO')
      MVLAF=INF(56)
      LA34F=INF(73)
      NSPSTA=INF(64)+1
      ALLOCATE(QTOT(1:NAF),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of QTOT')
C *** COMMON XYV,GVECT 
C *** CARTESIAN COMPONENTS OF LATTICE PARAMETERS
      CALL RREADF(IUNI,PARET,9)
C *** TRASFORMATION MATRIX CRYSTALLOGRAPHIC <=> PRIMITIVE CELL
      CALL RREADF(IUNI,W1R,9)
c *** INVERSE OF SYMMETRY OPERATOS
      CALL IREADF(IUNI,NINV,INF(2))
C *** MULTIPLICATION TABLE
      CALL IREADF(IUNI,MULTAB,48*48)
C ROTAtIONAL PART OF SYMMETRY OPERATORS IN CARTESIAN REFERENCE
      CALL RREADF(IUNI,XYV,9*INF(2))
C *** TRASLATIONAL PART OF SYMMETRY OPERATORS
      CALL RREADF(IUNI,TRASV,3*INF(2))
c *** print of CARTESIAN COMPONENTS OF LATTICE PARAMETERS
      IF(INF(10).NE.0)WRITE(IOUT,1130)((PARET(I,J),J=1,3),I=1,3)
1130  FORMAT(/' STRUCTURE AND SYMMETRY INFORMATION '//
     *' DIRECT LATTICE VECTOR COMPONENTS (BOHR)'//
     *T13,'X', T28,'Y',T43,'Z'/
     *' B1',3F15.9/' B2',3F15.9/' B3',3F15.9/)
C *** PRINT OF SYMMETRY OPERATORS
      MVF=INF(2)
      WRITE(IOUT,200)
      DO MV1=1,MVF,2
      MV9=MIN(MV1+1,MVF)
      WRITE(IOUT,201)(I,NINV(I),I=MV1,MV9)
      DO I=1,3
      WRITE(IOUT,202)((XYV(I,J,IV),J=1,3),TRASV(I,IV),IV=MV1,MV9)
      ENDDO
      ENDDO
C *** PRINT OF MULTIPLICATION TABLE
      WRITE(IOUT,100)
      CALL MATINT_D(MULTAB,MVF,48)
100   FORMAT(/' POINT GROUP MULTIPLICATION TABLE')
200   FORMAT(/' ***** SYMMOPS - TRANSLATORS IN BOHR')
201   FORMAT(/2('   NO.',I3,' INVERSE',I3,20X))
202   FORMAT(2(1X,3F8.3,2X,F8.3,5X))
      INF5=INF(5)+1
      INF79=INF(79)
      INF793=INF79*3
C *** INFORMATION ON DIRECT LATTICE - NUMBER OF DIRECT LATTICE VECTORS
      IF(INF79.GT.LIM006)CALL ERRNIC(0,INF79,ZNAMZ,
     *'TOO MANY G - INCREASE LIM006 TO')
C *** INFORMATION ON DIRECT LATTICE - NUMBER OF STARS OF DIRECT LATTICE VECTORS
      IF(INF(5).GT.LIM007)CALL ERRNIC(0,INF(5),ZNAMZ,
     *'TOO MANY G - INCREASE LIM007 TO')
C *** SQUARE OF MODULUS OF DIRECT LATTICE STARS OF VECTORS
      CALL RREADF(IUNI,GMODUS,INF5)
C *** CARTESIAN COORDINATES OF DIRECT LATTICE VECTORS
      CALL RREADF(IUNI,XG,INF793)
      CALL IREADF(IUNI,NM,INF5)
C *** NUMBER OF VETORS IN EACH STAR
      CALL IREADF(IUNI,MN,INF5)
C *** INVERSE VECTOR
      CALL IREADF(IUNI,NN1,INF79)
C *** COORDINATES OF LATTICE VECTORS IN CRYSTALLOGRAPHIC UNITS
      CALL IREADF(IUNI,LG,INF793)
C... PRINT THE STARS OF DIRECT LATTICE VECTORS
      NSTAR=INF(5)
c lprint(1)=3, number of stars of direct lattice vectors to be printed
      lprint(1)=3
      IF(INF(10).NE.0)THEN
      WRITE(IOUT,101)INF(79),INF(5),SQRT(GMODUS(NSTAR))
      JMAX=MIN(LPRINT(1),NSTAR)
      DO 14 LSH=1,JMAX
      N=MN(LSH)
      J=MN(LSH+1)
      M=J-N
      R=SQRT(GMODUS(LSH))
      WRITE(IOUT,102)LSH,M,R
      DO M=N,J-1
      WRITE(IOUT,103)M,NN1(M),(LG(I,M),I=1,3),(XG(I,M),I=1,3)
      ENDDO
14    CONTINUE
      ENDIF
101   FORMAT(/' NO.OF G VECTORS',I5,' STARS',I5,' RMAX',F11.5,' AU')
102   FORMAT(/' STAR N. ',I3,' N. OF VECTORS ',I3,' R=',1PE15.7/
     *'   LG I(LG) LX   LY   LZ         X           Y          Z')
103   FORMAT(5I5,4X,3F12.6,2I6)
204   FORMAT(/2X,A24)
c *** END PRINTING DIRECT LATTICE VECTORS
C *** COMMON BASATO
C *** NUCLEAR CHARGE
      CALL RREADF(IUNI,AZNUC,NAF)
C *** CARTESIAN CORDINATES OF ATOMS IN THE REFERENCE CELL
      CALL RREADF(IUNI,XA,NAF3)
C *** FORMAL CHARGE ATTRIBUTED TO EACH SHELL
      CALL RREADF(IUNI,CHE,LAF)
C *** ADJOINED GAUSSIAN OF EACH SHELL
      CALL RREADF(IUNI,EXAD,LAF)
C *** CARTESIAN COMPONENTS OF SHELL COORDINATES
      CALL RREADF(IUNI,XL,LAF3)
C *** PRIMITIVE GAUSSIAN EXPONENTS
      CALL RREADF(IUNI,EXX,NPRIM)
C *** CONTRACTION COEFFICIENTS OF GAUSSIANS IN CRYSTAL NORMALIZATION
      CALL RREADF(IUNI,C1,NPRIM)
      CALL RREADF(IUNI,C2,NPRIM)
      CALL RREADF(IUNI,C3,NPRIM)
      CALL RREADF(IUNI,CMAX,NPRIM)
C *** CONTRACTION COEFFICIENTS OF GAUSSIANS IN OLD NORMALIZATION 
      CALL RREADF(IUNI,C2W,NPRIM)
      CALL RREADF(IUNI,C3W,NPRIM)
C *** FORMAL ATOMIC NUMBER OF ATOMS (Z=MOD(NAT,100))
      CALL IREADF(IUNI,NAT,NAF)
C *** FIRST SHELL OF EACH ATOM
      CALL IREADF(IUNI,NSHPRI,NAF1)
C *** FIRST PRIMITIVE OF EACH SHELL
      CALL IREADF(IUNI,LAA,LAF1)
C *** NUMEBR OF PRIMITIVES IN EACH SHELL 
      CALL IREADF(IUNI,LAN,LAF)
C *** SHEL TYPE (0=S, 1=SP, 2=P, 3=D, 4=F)
      CALL IREADF(IUNI,LAT,LAF)
C *** NUMBER OF AO PER SHELL
      CALL IREADF(IUNI,LATAO,LAF)
C *** FIRST AO OF EACH SHELL
      CALL IREADF(IUNI,NDQ,LAF1)
C *** ATOMS TO WHICH EACH SHELL BELONGS
      CALL IREADF(IUNI,LATOAT,LAF)
      WRITE(IOUT,1072)
      DO NA=1,INF(24)
      IAM=NAT(NA)
      IA=MOD(IAM,100)
       WRITE(IOUT,1082)NA,IAM,SYMBAT(IA),(XA(K,NA),K=1,3)
      ENDDO
      WRITE(IOUT,1071)
1082  FORMAT(I4,2X,   I4,1X,A,1P,3E20.12)
1072  FORMAT(/' ATOM CARTESIAN COORDINATES (BOHR) - PRIMITIVE CELL'/
     *1X,79('*')/' *      ATOM',T22,'X',T42,'Y',T62,'Z'/1X,79('*'))
c *** PRINTING OF BASIS SET
      write(iout,1069)
1069  format(' VARIATIONAL BASIS SET'/)
      IPRAT(0:399)=.FALSE.
      WRITE(IOUT,1070)
      DO 430  NA=1,INF(24)
      IAM=NAT(NA)
      IA=MOD(IAM,100)
       WRITE(IOUT,1080)NA,SYMBAT(IA),(XA(K,NA),K=1,3)
      IF(IPRAT(IAM))GOTO 430
      IPRAT(IAM)=.TRUE.
      DO 420 I=NSHPRI(NA),NSHPRI(NA+1)-1
      LATI=LAT(I)+1
      K=NDQ(I)+1
      IF(LATI.EQ.1)THEN
      WRITE(IOUT,1090)K,NDN(LATI)
      ELSE
      WRITE(IOUT,1100)K,NDQ(I+1),NDN(LATI)
      ENDIF
      DO K=LAA(I),LAA(I+1)-1
      WRITE(IOUT,1120)EXX(K),C1(K),C2(K),C3(K)
      ENDDO
420   CONTINUE
430   CONTINUE
      WRITE(IOUT,1071)
1080  FORMAT(I4,1X,A,3F7.3)
1070  FORMAT(1X,79('*')/' LOCAL ATOMIC FUNCTIONS BASIS SET'/1X,79('*')/
     *'   ATOM  X(AU)  Y(AU)  Z(AU)    NO. TYPE  EXPONENT',
     *'  S COEF   P COEF   D/F/G COEF'/1X,79('*'))
1090  FORMAT(31X,I4,A)
1100  FORMAT(26X,I4,'-',I4,A)
1120  FORMAT(40X,1P,4E10.3)
1071  FORMAT(1X,79('*')/)
C *** COMMON INFPOT - INFORMATION ON ECP
      IF(INF(31).NE.0)THEN
      CALL IREADF(IUNI,IPSEUD,NAF)
      READ(IUNI,108)II,JJ,ITYPSE,(NPOT(I),I=1,II),(NBTYP(I),I=1,JJ),
     +    (NSOM(I),I=1,ITYPSE)
      READ(IUNI,109)(APOT(I),I=1,II),(CPOT(I),I=1,II)
C *** ECP PRINTING
      WRITE(IOUT,1111)
      JJ=0
      ITYPSE=0
      DO NA=1,NAF
      IF(IPSEUD(NA).EQ.0)CYCLE
      JJ=JJ+NANGECP
      ITYPSE=ITYPSE+1
      WRITE(IOUT,8730)MOD(NAT(NA),100),AZNUC(NA)
      IR=JJ-NANGECP
      K=NSOM(ITYPSE)
      DO J=1,NANGECP
      IR=IR+1
      NGAU=NBTYP(IR)
      IF(NGAU.EQ.0)CYCLE
      NB=MIN(NGAU,2)
      WRITE(IOUT,8750)XMOT(J),(APOT(I+K),CPOT(I+K),NPOT(I+K),I=1,NB)
      DO L=NB+1,NGAU,2
      LB=L+1
      IF(L.EQ.NGAU)LB=NGAU
      WRITE(IOUT,8760)(APOT(I+K),CPOT(I+K),NPOT(I+K),I=L,LB)
      ENDDO
      K=K+NGAU
      ENDDO
      ENDDO
      ENDIF
1111  FORMAT(/1X,79('*')/' *** PSEUDOPOTENTIAL INFORMATION ***'/
     *1X,79('*'))
8730  FORMAT(/' ATOMIC NUMBER',I4,', NUCLEAR CHARGE',F7.3/
     *' PSEUDOPOTENTIAL'/'   TYPE       EXPONENT      COEFF.    N',
     *'      EXPONENT      COEFF.    N')
 8750 FORMAT(2X,A6,4(F14.7,F13.7,I4))
 8760 FORMAT(8X,4(F14.7,F13.7,I4))
C *** COMMON MOLBAR - FOR MOLSPLIT ONLY
      IF(INF(34).NE.0)THEN
      READ(IUNI,108)(N1MOL(I),I=1,INF(92)+1)
      CALL RREADF(IUNI,XBAR,INF(92)*3)
      WRITE(IOUT,110)INF(92)
      WRITE(IOUT,111)(J,(XBAR(I,J),I=1,3),J=1,INF(92))
      WRITE(IOUT,112)(I,N1MOL(I),I=1,INF(92))
111   FORMAT(/' COORDINATES OF THE BARYCENTRES OF THE MOLECULES'//
     *(I6,3E20.12))
110   FORMAT(/' LATTICE OF ',I3,' MOLECULES NOT INTERACTING')
112   FORMAT(/' FIRST ATOM OF THE MOLECULE:'/(8(I5,I4,',')))
      ENDIF
C *** ELECTRONIC CHARGE ACCORDING TO MULLIKEN POPULATION ANALYSIS
      CALL RREADF(IUNI,QTOT,NAF)
      WRITE(IOUT,1079)
1079  FORMAT(/' MULLIKEN ATOMIC CHARGES'/)
1081  FORMAT(2(I4,1X,A,' ZNUC',F7.3,' QTOT',F7.3))
      WRITE(IOUT,1081)(NA,SYMBAT(MOD(NAT(NA),100)),AZNUC(NA),QTOT(NA),
     *NA=1,NAF)
106   FORMAT(8I10)
107   FORMAT(80A1)
203   FORMAT(/13(' CRY06')//' INFORMATION FROM CRYSTAL06'//1X,80A1//
     *' TYPE OF CALCULATION : ',A24/)
108   FORMAT(/(8I10))
109   FORMAT(/(4E20.13))
C *** COMMON LAVLAF
      MVLAF1=MVLAF+1
      LAF1=LAF+1
c     CALL IREADF(IUNI,LA34F,1)
      ALLOCATE(ncf(1:mvlaf1),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of ncf')
      ALLOCATE(inzvlb(1:mvlaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of inzvlb')
      ALLOCATE(idime(1:mvlaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of idime')
      ALLOCATE(idimf(1:mvlaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of idimf')
      ALLOCATE(idmcou(1:mvlaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of idmcou')
C *** SYMMETRY INFORMATION - FOR PROGRAMMERS ONLY
      CALL IREADF(IUNI,NCF,MVLAF1)
      INF793=NCF(MVLAF1)
      CALL IREADF(IUNI,INZVLB,MVLAF)
      CALL IREADF(IUNI,IDIME,MVLAF)
      CALL IREADF(IUNI,IDIMF,MVLAF)
      CALL IREADF(IUNI,IDMCOU,MVLAF)
      ALLOCATE(jncdu(1:mvlaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of jncdu')
      ALLOCATE(irof(1:mvlaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of irof')
      ALLOCATE(jrof(1:mvlaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of jrof')
      ALLOCATE(lav(1:mvf*laf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of lav')
      ALLOCATE(mgnav(1:mvf*naf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of mgnav')
      ALLOCATE(la3(1:inf793),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of la3')
      ALLOCATE(la4(1:inf793),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of la4')
      CALL IREADF(IUNI,LAV,MVF*LAF)
      CALL IREADF(IUNI,MGNAV,MVF*NAF)
      CALL IREADF(IUNI,LA3,INF793)
      CALL IREADF(IUNI,LA4,INF793)
      CALL IREADF(IUNI,JNCDU,MVLAF)
      CALL IREADF(IUNI,IROF,MVLAF)
      CALL IREADF(IUNI,JROF,MVLAF)
      READ(IUNI,*)isizeiccs3,isizeiccs1,isizeilana
      ALLOCATE(iccat(1:naf1),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of iccat')
      ALLOCATE(iccs3(1:isizeiccs3),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of iccs3')
      ALLOCATE(kvrsp(1:laf*laf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of kvrsp')
      ALLOCATE(iccv(1:laf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of iccv')
      ALLOCATE(icct(1:laf1),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of icct')
      ALLOCATE(icc(1:laf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,IERROR,ZNAMZ,'allocation of icc')
      CALL IREADF(IUNI,ICCAT,NAF1)
      CALL IREADF(IUNI,ICCS3,isizeiccs3)
      CALL IREADF(IUNI,KVRSP,LAF*LAF)
      CALL IREADF(IUNI,iccv,LAF)
      ICCLAF=ICCAT(NAF1)
      ALLOCATE(ina12(1:icclaf-1),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf-1,ZNAMZ,'allocation of ina12')
      ALLOCATE(nlana(1:icclaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf,ZNAMZ,'allocation of nlana')
      CALL IREADF(IUNI,INA12,ICCLAF-1)
      CALL IREADF(IUNI,NLANA,ICCLAF)
      CALL IREADF(IUNI,ICC,LAF)
      CALL IREADF(IUNI,ICCT,LAF1)
      ICCLAF=ICCT(LAF1)
      ALLOCATE(ila12t(1:icclaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf,ZNAMZ,'allocation of ila12t')
      ALLOCATE(jpoint(1:icclaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf,ZNAMZ,'allocation of jpoint')
      ALLOCATE(iccs1(1:iSIZEICCS1),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,iSIZEICCS1,ZNAMZ,
     *'allocation of iccs1')
      ALLOCATE(nla21t(1:icclaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf,ZNAMZ,'allocation of nla21t')
      ALLOCATE(nnnc(1:icclaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf,ZNAMZ,'allocation of nnnc')
      ALLOCATE(nnnc2(1:icclaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf,ZNAMZ,'allocation of nnnc2')
      ALLOCATE(la34x(1:icclaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf,ZNAMZ,'allocation of la34x')
      ALLOCATE(la34v(1:icclaf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,icclaf,ZNAMZ,'allocation of la34v')
      ALLOCATE(ilana(1:iSIZEILANA),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,iSIZEILANA,ZNAMZ,
     *'allocation of ilana')
      CALL IREADF(IUNI,ILA12T,ICCLAF)
      CALL IREADF(IUNI,JPOINT,ICCLAF)
      CALL IREADF(IUNI,ICCS1,iSIZEICCS1)
      CALL IREADF(IUNI,NLA21T,ICCLAF)
      CALL IREADF(IUNI,NNNC,ICCLAF)
      CALL IREADF(IUNI,NNNC2,ICCLAF)
      CALL IREADF(IUNI,LA34X,ICCLAF)
      CALL IREADF(IUNI,LA34V,ICCLAF)
      CALL IREADF(IUNI,ILANA,iSIZEILANA)
!     read(iuni,106)iSIZENGSHG,iSIZENQGSHG,iSIZENNGI,
!    *isizenshgi!,iSIZEMGV,iSIZEINO,iSIZEINOIV
      NNGIDMF=INF(39)*INF(133)
      ALLOCATE(nshg(1:la34f),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,la34f,ZNAMZ,'allocation of nshg')
      ALLOCATE(ngshg(1:la34f*inf(37)),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,la34f*inf(37),ZNAMZ,
     *'allocation of ngshg')
      ALLOCATE(nqgshg(1:la34f*inf(145)),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,la34f*inf(145),ZNAMZ,
     *'allocation of nqgshg')
      ALLOCATE(nstatg(1:mvlaf1),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,mvlaf1,ZNAMZ,'allocation of nqgshg')
      ALLOCATE(nngi(1:nngidmf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,nngidmf,ZNAMZ,'allocation of nngi')
      ALLOCATE(nshgi(1:nngidmf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,nngidmf,ZNAMZ,'allocation of nshgi')
      CALL IREADF(IUNI,NSHG,LA34F)
      CALL IREADF(IUNI,NGSHG,la34f*inf(37))
      CALL IREADF(IUNI,NQGSHG,la34f*inf(145))
C *** COMMON INOETC
      CALL IREADF(IUNI,NSTATG,MVLAF1)
      CALL IREADF(IUNI,NNGI,NNGIDMF)
      CALL IREADF(IUNI,NSHGI,NNGIDMF)
c rotmatrix
!     LIMROT1=INF(183)
!     ALLOCATE(TTO(LIMROT1),STAT=IERR)
!     IF(IERR.NE.0)CALL ERRNIC(0,LIMROT1,ZNAMZ,'TTO ALLOCATION')
!     ALLOCATE(MMO(LIMROT1),STAT=IERR)
!     IF(IERR.NE.0)CALL ERRNIC(0,LIMROT1,ZNAMZ,'MMO ALLOCATION')
!     ALLOCATE(MMOM(LIMROT1),STAT=IERR)
!     IF(IERR.NE.0)CALL ERRNIC(0,LIMROT1,ZNAMZ,'MMOM ALLOCATION')
      NPGT=INF(12)
      NFGT=INF(11)
      ALLOCATE(sg(1:npgt),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,npgt,ZNAMZ,'allocation of sg')
      ALLOCATE(fg(1:nfgt*NSPSTA),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,nfgt,ZNAMZ,'allocation of fg')
      ALLOCATE(pg(1:npgt*NSPSTA),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,npgt,ZNAMZ,'allocation of pg')
      ALLOCATE(a(1:ndf*ndf),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,ndf*ndf,ZNAMZ,'allocation of a')
c iky computed to print hamiltonian matrix
      ALLOCATE(iky(1:ndf+1),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,ndf,ZNAMZ,'allocation of iky')
c fort.3  overlap matrix
cro      CALL RREADF(IUNI,SG,NPGT)
      CALL RREADF(IUNI,SG,NFGT)
c *** UNRESTRICTED OPEN SHELL:
C *** FOCK MATRIX: 2 ARRAYS, ALPHA and BETA
c fort.11 Fock/KS matrix
      CALL RREADF(IUNI,FG,NFGT)
      IF(NSPSTA.NE.1)CALL RREADF(IUNI,FG(NFGT+1),NFGT)
c *** UNRESTRICTED OPEN SHELL:
C *** DENSITY MATRIX: 1 ARRAY, ALPHA+BETA, ALFA-BETA (total, spin)
c fort.13 density matrix
      CALL RREADF(IUNI,PG,NPGT*NSPSTA)
C
c *** printing of density matrix
C
      write(iout,180)
180   format(/' DENSITY MATRIX DIRECT LATTICE - G=0'/)
c
      MGMG=1
c matpmn - density matrix direct space - mgmg direct lattice vector
c modify the value of MGMG to print the density matrix relative to the
c first mgmg direct lattice vectors
      NDF=INF(7)
      LAF=INF(20)
      NSTAP=0
      DO ISIGMA=1,NSPSTA
      WRITE(IOUT,204)NSPIN(ISIGMA)
      A(1:NDF*NDF)=0._FLOAT
      K=0
      DO 124 L1=1,LAF
      LA1=LATAO(L1)
      DO 123 LL2=ICCT(L1)+1,ICCT(L1+1)
      L2=ILA12T(LL2)
      KSH=IDIME(JPOINT(LL2))
      IF(KSH.EQ.0)GOTO 123
      IPPOP=LA34V(LL2)
      INZ=NGSHG(IPPOP+1)+1
      DO 1244 MG=INZ,NGSHG(IPPOP+KSH+1)
      IF(NQGSHG(MG).NE.MGMG)GOTO 1244
      LA2=LATAO(L2)
      II=NDQ(L2)+K
      M=(MG-INZ)*LA1*LA2+NNNC(LL2)
      DO I=1,LA1
      DO J=1,LA2
      M=M+1
      A(II+J)=PG(M+NSTAP)
      ENDDO
      II=II+NDF
      ENDDO
      GOTO 123
 1244 CONTINUE
123   CONTINUE
      K=LA1*NDF+K
124   CONTINUE
      DO M=1,NDF,10
      K=MIN(M+9,NDF)
      WRITE(IOUT,40) (J,J=M,K)
      K=K-M
      INZ1=M
      DO I=1,NDF
      INZ2=INZ1+K
      WRITE(IOUT,50)I,(A(J),J=INZ1,INZ2)
      INZ1=INZ1+NDF
      ENDDO
      ENDDO
      NSTAP=NPGT
      ENDDO
50    FORMAT(I4,3X,1P,10E12.4)
40    FORMAT(/3X,10(8X,I4)/)
      MGMG=1
C
c matfg ini - hamiltonian matrix in direct space - mgmg direct lattice vector
c modify the value of MGMG to print the density matrix relative to the
c first mgmg direct lattice vectors
C
      write(iout,181)
181   format(/' HAMILTONIAN MATRIX DIRECT LATTICE - G=0')
      NDF=INF(7)
      LAF=INF(20)
      NSTAP=0
      DO ISIGMA=1,NSPSTA
      IF(NSPSTA.NE.1)WRITE(IOUT,204)MSPIN(ISIGMA)
      II=0
      DO I=1,NDF+1
      IKY(I)=II
      II=II+I
      ENDDO
      A(1:IKY(NDF+1))=0._FLOAT
      DO L1=1,LAF
      LA1=LATAO(L1)
      INZ1=NDQ(L1)+1
      IFN1=NDQ(L1+1)
      DO L2=L1,LAF
      LL2=KVRSP(ICCV(L1)+L2)
      IF(LL2.EQ.0)CYCLE
      KSH=IDIMF(JPOINT(LL2))
      IF(KSH.EQ.0)CYCLE
      IPPOP=ICCS3(LL2)
      MSTART=NNNC2(LL2)
      INZ=NGSHG(IPPOP+1)+1
      IF(L1.EQ.L2)THEN
      MMULT=IKY(LA1+1)
      ELSE
      LA2=LATAO(L2)
      INZ2=NDQ(L2)
      MMULT=LA1*LA2
      ENDIF
      DO MG=INZ,NGSHG(IPPOP+KSH+1)
      IF(NQGSHG(MG).NE.MGMG)CYCLE
      M=(MG-INZ)*MMULT+MSTART
      IF(L1.EQ.L2)THEN
      DO I=INZ1,IFN1
      II=IKY(I)
      DO J=INZ1,I
      M=M+1
      A(II+J)=FG(M)
      ENDDO
      ENDDO
      ELSE
      DO I=1,LA2
      II=IKY(I+INZ2)
      MM=M+I
      DO J=INZ1,IFN1
      A(II+J)=FG(MM+NSTAP)
      MM=MM+LA2
      ENDDO
      ENDDO
      ENDIF
      EXIT
      ENDDO
      ENDDO
      ENDDO
      II=1
      DO M=1,NDF,10
      I=MIN(M+9,NDF)
      WRITE(IOUT,40)(J,J=M,I)
      DO I=M,NDF
      INZ1=IKY(I)+II
      INZ2=INZ1+MIN(I-M,9)
      WRITE(IOUT,50)I,(A(J),J=INZ1,INZ2)
      ENDDO
      II=II+10
      ENDDO
      NSTAP=NFGT
      ENDDO
c matfg end
c
      if(inf(139).ne.0)THEN 
      WRITE(IOUT,1068)
1068  FORMAT(/' LOCALIZATION DATA')
      DO ISIGMA=1,NSPSTA
      IF(NSPSTA.NE.1)WRITE(IOUT,204)MSPIN(ISIGMA)
      CALL READLOC(NDF,NAF,IUNI)
      ENDDO
      ENDIF
      CLOSE(IUNI)
      DEALLOCATE(iky)
      DEALLOCATE(a)
      DEALLOCATE(pg)
      DEALLOCATE(fg)
      DEALLOCATE(sg)
      DEALLOCATE(nshgi)
      DEALLOCATE(nngi)
      DEALLOCATE(nstatg)
      DEALLOCATE(nqgshg)
      DEALLOCATE(ngshg)
      DEALLOCATE(nshg)
      DEALLOCATE(ilana)
      DEALLOCATE(la34v)
      DEALLOCATE(la34x)
      DEALLOCATE(nnnc2)
      DEALLOCATE(nnnc)
      DEALLOCATE(nla21t)
      DEALLOCATE(iccs1)
      DEALLOCATE(jpoint)
      DEALLOCATE(ila12t)
      DEALLOCATE(nlana)
      DEALLOCATE(ina12)
      DEALLOCATE(icc)
      DEALLOCATE(icct)
      DEALLOCATE(iccv)
      DEALLOCATE(kvrsp)
      DEALLOCATE(iccs3)
      DEALLOCATE(iccat)
c     seg fault if executed
      DEALLOCATE(la4)
      DEALLOCATE(la3)
      DEALLOCATE(mgnav)
      DEALLOCATE(lav)
      DEALLOCATE(jrof)
      DEALLOCATE(irof)
      DEALLOCATE(jncdu)
      DEALLOCATE(idmcou)
      DEALLOCATE(idimf)
      DEALLOCATE(idime)
      DEALLOCATE(inzvlb)
c     seg fault if executed
      DEALLOCATE(ncf)
      DEALLOCATE(QTOT)
c     DEALLOCATE(PAR)
c     DEALLOCATE(ITOL)
c     DEALLOCATE(INF)
c     REWIND IUNI
      RETURN
      END
      SUBROUTINE READLOC(NDF,NAF,IU)
      Use numbers
      Use parame_module
      Use parinf_module
      Use gvect_module
      Implicit real( float ) ( a-h, o-z )
      REAL(FLOAT),DIMENSION(:,:,:),ALLOCATABLE :: wfun
      REAL(FLOAT),DIMENSION(:,:),ALLOCATABLE :: wcent
      REAL(FLOAT),DIMENSION(:),ALLOCATABLE :: poploc
      INTEGER,DIMENSION(:),ALLOCATABLE :: lbands
      CHARACTER(LEN=7) :: ZNAMZ='READLOC'
      READ(IU,106)MMGG,MMMG0,LIMCLL,NOCC
C +++ MMGG, MMMG0 used by external program
      ALLOCATE(LBANDS(1:NOCC),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,NOCC,ZNAMZ,
     *'MEMORY ALLOCATION FOR LBANDS - FAILED')
      CALL IREADF(IU,LBANDS,NOCC)
      WRITE(IOUT,1005)(LBANDS(I),I=1,NOCC)
C +++ READ WANNIER FUNCTIONS 
      ALLOCATE(WFUN(1:NDF,1:NOCC,1:MMMG0),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,MMMG0*NOCC*NDF,ZNAMZ,
     *'MEMORY ALLOCATION FOR WFUN - FAILED')
      CALL RREADF(IU,WFUN,MMMG0*NOCC*NDF)
C +++ print WF for the first 3 lattice vectors
      DO MG=1,MIN(3,MMMG0)
      WRITE(IOUT,1000)(LG(I,MG),I=1,3)
      CALL CMACOL(WFUN(1,1,MG),NDF,1,NOCC,1)
      ENDDO
C +++ READ WFS CENTROIDS 
      ALLOCATE(WCENT(1:3,1:NOCC),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,NOCC*3,ZNAMZ,
     *'MEMORY ALLOCATION FOR WCENT - FAILED')
      CALL RREADF(IU,WCENT,NOCC*3)
      WRITE(IOUT,1004)
      DO IC=1,NOCC
      WRITE(IOUT,1006)IC,(WCENT(N,IC),N=1,3)
      ENDDO
1004  FORMAT(//' WANNIER',T26,'CENTROID''S COORDINATES (BOHR):'/)
1006  FORMAT(I6,3E20.12)
C +++ READ ATOMIC POPULATION FROM UNIT 80
      LNM=NAF*NOCC*MMMG0
      ALLOCATE(POPLOC(LNM),STAT=IERROR)
      IF(IERROR.NE.0)CALL ERRNIC(0,LNM,ZNAMZ,
     *'MEMORY ALLOCATION FOR POPLOC - FAILED')
      CALL RREADF(IU,POPLOC,LNM)
      DEALLOCATE(POPLOC)
      DEALLOCATE(WCENT)
      DEALLOCATE(WFUN)
      DEALLOCATE(LBANDS)
      RETURN
1000  FORMAT(/'G = (',3I4,')')
 1005 FORMAT(/' WANNIER FUNCTIONS - LIST OF ACTIVE BANDS'/
     *(10I7/))
106   FORMAT(8I10)
      END
      SUBROUTINE CMACOL(A,N,NINIZ,NIFIN,LATV)
      Use numbers
      Use parame_module
      Use parinf_module
      Implicit real( float ) ( a-h, o-z )
      DIMENSION A(*),FIT(10)
      N10=N*10
      IF(LATV.EQ.0)THEN
      N2=N+N
      MBASE=(NINIZ-1)*N2
      DO M=NINIZ,NIFIN,5
      K=MIN(M+4,NIFIN)
      WRITE(IOUT,40)(J,J,J=M,K)
      DO I=1,N
      NBASE=MBASE+I+I
      L=0
      DO L1=M,K
      FIT1=A(NBASE-1)
      FIT2=A(NBASE)
      IF(ABS(FIT1).LT.(1.E-7))FIT1=0.
      IF(ABS(FIT2).LT.(1.E-7))FIT2=0.
      FIT(L+1)=FIT1
      FIT(L+2)=FIT2
      NBASE=NBASE+N2
      L=L+2
      ENDDO
      WRITE(IOUT,50)I,(FIT(L1),L1=1,L)
      ENDDO
      MBASE=MBASE+N10
      ENDDO
      ELSE
      MBASE=(NINIZ-1)*N
      DO M=NINIZ,NIFIN,10
      K=MIN(M+9,NIFIN)
      WRITE(IOUT,40)(J,J=M,K)
      DO I=1,N
      NBASE=MBASE+I
      L=0
      DO L1=M,K
      L=L+1
      FIT1=A(NBASE)
      IF(ABS(FIT1).LT.(1.E-7))FIT1=0.
      FIT(L)=FIT1
      NBASE=NBASE+N
      ENDDO
      WRITE(IOUT,50)I,(FIT(L1),L1=1,L)
      ENDDO
      MBASE=MBASE+N10
      ENDDO
      ENDIF
      RETURN
50    FORMAT(I5,3X,10G12.4)
40    FORMAT(/4X,10(8X,I4)/)
      END
C
      SUBROUTINE KREDIN
C *** THIS ROUTINE WAS DERIVED FROM BMAT/SMAT ROUTINE
C *** READS EIGENVECTORS IN IBZ AND GENERATE EIGENVECTORS
C *** IN THE FULL BRILLOUIN ZONE, APPLYING SYMMETRY
C *** OPERATORS AND TIME REVERSAL SYMMETRY
      USE NUMBERS
      Use parinf_module
      Use retic_module
      Use epesi_module
      Use xyvdim_module
      Use text_module
      IMPLICIT REAL(FLOAT) (A-H,O-Z)
      LOGICAL exist
!     Real( float ),Dimension(:),Allocatable :: a
      Real( float ),Dimension(:),Allocatable :: ar
      Real( float ),Dimension(3):: akxyz
      Integer,Dimension(:),Allocatable :: norder
      Integer,Dimension(3):: kr
      CHARACTER(LEN=6) :: ZNAMZ='kredin'
      MVF=INF(2)
      NDF=INF(7)
      NSPSTA=INF(64)+1
      ndfsqr2=ndf*ndf*2
      Allocate(ar(ndfsqr2),Stat=ierror )
      If(ierror.ne.0)call errnic(0,ndfsqr2,znamz,'ar Memory allocation')
      NBAND=NDF*NSPSTA
      IO30=IUNIT(30)
!
      INQUIRE(FILE='KRED.DAT',EXIST=EXIST)
      IF(.NOT.EXIST)THEN
         CALL ERRNIC(0,NDF,ZNAMZ,
     *'FILE KRED.DAT NOT FOUND - EIGENVECTORS')
         ELSE
           OPEN(UNIT=IO30,FILE='KRED.DAT',FORM='FORMATTED',
     *STATUS='OLD')
         ENDIF
      READ(IO30,101)IS1,IS2,IS3,NKF,BRET
      IF(NKF.GT.LIM001)CALL ERRNIC(0,NKF,ZNAMZ,
     *'TOO MANY K POINTS - INCREASE LIM001 TO')
      IS1IS2=IS1*IS2
      Allocate(norder(is1*is2*is3),Stat=ierror )
      If(ierror.ne.0)call errnic(0,is1*is2*is3,znamz,
     *'norder Memory allocation')   
C NKF          NUMBER OF K POINTS IN IBZ
C INF(1)       NUMBER OF SYMMETRY OPERATORS EXTENDED TO INVERSION
C INF(2)       NUMBER OF SYMMETRY OPERATORS
C INF(7)       NUMBER OF BASIS FUNCTIONS (AO)
C INF(10)      TRANSLATIONAL SYMMETRY: 0 (0D) 1 (1D) 2 (2D) 3 (3D)
C IS1(IS2,IS3) MONKHORST NET SHRINKING FACTORS 
C INF(64)      0 (RESTRICTED CLOSED SHELL) / 1 (UNRESTRICTED OPEN SHELL)
      NTOT=(INF(64)+1)*INF(7)*NKF
C *** COORDINATES OF K POINTS IN LATTICE VECTORS UNITS
      CALL IREADF(IO30,JJ,3*NKF)
C *** TYPE OF K POINT: Complex (0), REAL (1)
      CALL IREADF(IO30,LATVRS,NKF)
C *** SYMMMETRY OPERATORS IN lATTICE VECTORS UNITS
      CALL IREADF(IO30,IRR,9*48)
C *** GEOMETRICAL WEIGHT OF K POINTS 
      CALL RREADF(IO30,WPJ,NKF)
      IF(NTOT.GT.LIM018)CALL ERRNIC(0,NTOT,ZNAMZ,
     *'N. OF K POINTS x BANDS TOO LARGE - INCREASE LIM018 TO')
C EIGENVALUES - UNREStRICTED: ALPHA(ALL K POINTS),BETA(ALL K POINTS)
      CALL RREADF(IO30,ENE,NTOT)
C ALFA - WEIGHT OF THE EIGENVALUES - COMPUTED FROM FERMI ENERGY CALCULATION
      CALL RREADF(IO30,ALFA,NTOT)
      WRITE(IOUT,100)
      IF(INF(10).NE.0)THEN
      WRITE(IOUT,1121)IS1,IS2,IS3,NKF,IS1,
     * (NK,TIPO(LATVRS(NK)),(JJ(I,NK),I=1,3),NK=1,NKF)
1121  FORMAT(' SHRINKING FACTOR (MONKHORST NET)',T35,3I3,
     *T46,'NUMBER OF K POINTS IN THE IBZ',T76,I4/
     *1X,79('*')/' *** K POINTS COORDINATES (OBLIQUE COORDINATES IN UNIT
     *S OF IS =',I3,')'/(4(I4,'-',A1,'(',3I3,')')))
      WRITE(IOUT,1103)(WPJ(NK),NK=1,NKF)
1103  FORMAT(/' GEOMETRICAL WEIGHT OF K POINTS - MONKHORST'//(08F10.5))
1104  FORMAT(/' WEIGHT OF K POINTS FOR EACH BAND'/)
1105  FORMAT(' BAND N',I4/(15F7.4))
      NBAND=NSPSTA*INF(7)
      NDINZ=NKF*NBAND
      WRITE(IOUT,1104)
      DO N=1,NBAND
      WRITE(IOUT,1105)N,(ALFA(NKN),NKN=N,NDINZ,NBAND)
      ENDDO
      WRITE(IOUT,1129)
      DO I=1,3
      WRITE(IOUT,1130)(BRET(I,J),J=1,3)
      ENDDO
      ENDIF
      JNIZ=0
      NBAND=NDF*NSPSTA
      DO ISIGMA=1,NSPSTA
      IF(NSPSTA.NE.1)WRITE(IOUT,204)MSPIN(ISIGMA)
      INIZ=JNIZ
      DO  K=1,NKF
      LATW=LATVRS(K)
      WRITE(IOUT,105)K,JJ(1,K),JJ(2,K),JJ(3,K),(ENE(INIZ+I),
     *I=1,NDF)
      INIZ=INIZ+NBAND
      ENDDO
      JNIZ=NDF
      ENDDO
105   FORMAT(/' EIGENVALUES (AU) - K=',I4,' (',3I3,')',1P/(10E12.4))
1129  FORMAT(/' RECIPROCAL LATTICE VECTORS COMPONENTS. (A.U.) '//
     *6X,1(2X,'X',12X,'Y',12X,'Z',12X))
1130  FORMAT(3F13.7,2X,3F13.7)
100   FORMAT(39(' K')//' EIGENVALUES AND RELATED INFORMATION ARE READ',
     *' FROM FILE KRED.DAT (FULL BZ)'//39(' K')/)
101   FORMAT(4I4/1P,(3E21.13))
!
      IS10=IS1*16
      IS20=IS2*16
      IS30=IS3*16
      VRSISS=1._FLOAT/IS1
      INIZ=0
      IF(INF(10).EQ.0)THEN
      WRITE(IOUT,205)
      ELSE
      WRITE(IOUT,104)
      WRITE (IOUT, 200)
      DO MV1=1,MVF,4
      MV9=MIN(MV1+3,MVF)
      WRITE(IOUT,201)(I,I=MV1,MV9)
      DO I=1,3
      WRITE(IOUT,202)((IRR(I,J,IV),J=1,3),IV=MV1,MV9)
      ENDDO
      ENDDO
      ENDIF
200   FORMAT(/' SYMMOPS - RECIPROCAL LATTICE'/)
201   FORMAT(/4(7X,' NO.',I3,6X))
202   FORMAT(4(3X,3I4,5X))
204   FORMAT(/2X,A24)
      JNIZ=0
      DO ISIGMA=1,NSPSTA
      NORDER(1:IS1IS2*IS3)=0
      IF(NSPSTA.GT.1)WRITE(IOUT,204)MSPIN(ISIGMA)
      INIZ=JNIZ
      DO 1010 K=1,NKF
      LATW=LATVRS(K)
      WRITE(IOUT,115)K,TIPO(LATW),JJ(1,K),JJ(2,K),JJ(3,K),(ENE(INIZ+I),
     *I=1,NDF)
104   FORMAT(//' HAMILTONIAN EIGENVECTORS IN THE FULL BRILLOUIN ZONE'/)
205   FORMAT(//' HAMILTONIAN EIGENVECTORS'/) 
115   FORMAT(/' EIGENVALUES - K=',I4,'-',A1,'(',3I3,') (IBZ)',1P/
     *(10E12.4))
106   FORMAT(/' EIGENVECTORS - K=',I4,'-',A1,'(',3I3,') - SYMMOP ',I4,
     *' (',3I3,') (COORDINATES UNITS IS=',I3,')'/
     *'(',3f12.8,' A.U.)')
107   FORMAT(/' TIME REVERSAL SYMMETRY APPLIED ')
      INIZ=INIZ+NBAND
      JR1=JJ(1,K)
      JR2=JJ(2,K)
      JR3=JJ(3,K)
      IF(LATVRS(K).NE.0)GOTO 889
C........................ COMPLEX K-POINT ............................
      DO 776 MV=1,MVF
      MR1=MOD(IRR(1,1,MV)*JR1+IRR(1,2,MV)*JR2+IRR(1,3,MV)*JR3+IS10,IS1)
      MR2=MOD(IRR(2,1,MV)*JR1+IRR(2,2,MV)*JR2+IRR(2,3,MV)*JR3+IS20,IS2)
      MR3=MOD(IRR(3,1,MV)*JR1+IRR(3,2,MV)*JR2+IRR(3,3,MV)*JR3+IS30,IS3)
      KR(1)=MR1
      KR(2)=MR2
      KR(3)=MR3
      NREC=MR1+1+MR2*IS1+MR3*IS1IS2
      IF(NORDER(NREC).NE.0)GOTO 776
      DO I=1,3
      AKXYZ(I)=(BRET(1,I)*KR(1)+BRET(2,I)*KR(2)+BRET(3,I)*KR(3))*VRSISS
      ENDDO
      WRITE(IOUT,106)K,TIPO(LATW),JR1,JR2,JR3,MV,KR,IS1,AKXYZ
!     CALL IREADF(IO30,IRR(1,1,MV),9)
      CALL IREADF(IO30,KR,3)
Check: k generated vs k written by CRYSTAL03
      NORDER(NREC)=2
      CALL RREADF(IO30,AR,NDF*NDF*(2-LATW))
      CALL CMACOL(AR,NDF,1,NDF,LATW)
776   CONTINUE
C *** NO TIME REVERSAL SYMMETRY - SKIP
      IF(MVF.EQ.INF(1))GOTO 777
C *** TIME REVERSAL SYMMETRY APPLIED
      JR1=MOD(IS1-JR1,IS1)
      JR2=MOD(IS2-JR2,IS2)
      JR3=MOD(IS3-JR3,IS3)
      DO 778 MV=1,MVF
      MR1=MOD(IRR(1,1,MV)*JR1+IRR(1,2,MV)*JR2+IRR(1,3,MV)*JR3+IS10,IS1)
      MR2=MOD(IRR(2,1,MV)*JR1+IRR(2,2,MV)*JR2+IRR(2,3,MV)*JR3+IS20,IS2)
      MR3=MOD(IRR(3,1,MV)*JR1+IRR(3,2,MV)*JR2+IRR(3,3,MV)*JR3+IS30,IS3)
      KR(1)=MR1
      KR(2)=MR2
      KR(3)=MR3
      NREC=MR1+1+MR2*IS1+MR3*IS1IS2
      IF(NORDER(NREC).NE.0)GOTO 778
      DO I=1,3
      AKXYZ(I)=(BRET(1,I)*KR(1)+BRET(2,I)*KR(2)+BRET(3,I)*KR(3))*VRSISS
      ENDDO
      WRITE(IOUT,107)
      WRITE(IOUT,106)K,TIPO(LATW),JJ(1,K),JJ(2,K),JJ(3,K),MV,KR,
     *IS1,AKXYZ
!     CALL IREADF(IO30,IRR(1,1,MV),9)
      CALL IREADF(IO30,KR,3)
      NORDER(NREC)=2
      CALL RREADF(IO30,AR,NDF*NDF*(2-LATW))
      CALL CMACOL(AR,NDF,1,NDF,LATW)
778   CONTINUE
777   GOTO 1010
C.......................... REAL K-POINT .............................
889   CONTINUE
      DO 886 MV=1,MVF
      MR1=MOD(IRR(1,1,MV)*JR1+IRR(1,2,MV)*JR2+IRR(1,3,MV)*JR3+IS10,IS1)
      MR2=MOD(IRR(2,1,MV)*JR1+IRR(2,2,MV)*JR2+IRR(2,3,MV)*JR3+IS20,IS2)
      MR3=MOD(IRR(3,1,MV)*JR1+IRR(3,2,MV)*JR2+IRR(3,3,MV)*JR3+IS30,IS3)
      KR(1)=MR1
      KR(2)=MR2
      KR(3)=MR3
      NREC=MR1+1+MR2*IS1+MR3*IS1IS2
      IF(NORDER(NREC).NE.0)GOTO 886
      DO I=1,3
      AKXYZ(I)=(BRET(1,I)*KR(1)+BRET(2,I)*KR(2)+BRET(3,I)*KR(3))*VRSISS
      ENDDO
      WRITE(IOUT,106)K,TIPO(LATW),JR1,JR2,JR3,MV,KR,IS1,AKXYZ
      CALL IREADF(IO30,KR,3)
      print 1943,'kr',kr
1943  format(a,9i3)
      NORDER(NREC)=1
      CALL RREADF(IO30,AR,NDF*NDF*(2-LATW))
      CALL CMACOL(AR,NDF,1,NDF,LATW)
886   CONTINUE
1010  CONTINUE
      INIZ=NDF
      ENDDO
      CLOSE(IO30)
      RETURN
      END
      SUBROUTINE IREADF(JOUT,IA,N)
      USE NUMBERS
      IMPLICIT REAL(FLOAT) (A-H,O-Z)
      DIMENSION IA(N)
      READ(JOUT,*)IA
      RETURN
      END 
      SUBROUTINE RREADF(JIN,A,N)
      USE NUMBERS
      IMPLICIT REAL(FLOAT) (A-H,O-Z)
      DIMENSION A(N)
      READ(JIN,*)A
      RETURN
      END
      SUBROUTINE ERRNIC(IERR,INIC,NAMZ,MZSS)
      USE NUMBERS
      USE PARINF_MODULE
      IMPLICIT REAL(FLOAT) (A-H,O-Z)
      CHARACTER*(*) NAMZ
      CHARACTER*(*) MZSS 
      IF(IERR.EQ.0)THEN
       WRITE(IOUT,1)TRIM(NAMZ),TRIM(MZSS),INIC
1      FORMAT(' ERROR **** ',A,' **** ',A,I10)
       STOP
      ELSE IF(IERR.EQ.1) THEN
       WRITE(IOUT,2)TRIM(NAMZ),TRIM(MZSS),INIC
2      FORMAT(' WARNING **** ',A,' **** ',A,I10)
      ELSE
       WRITE(IOUT,3)TRIM(NAMZ),TRIM(MZSS),INIC
3      FORMAT(' INFORMATION **** ',A,' **** ',A,I10)
      ENDIF
      RETURN
      END
      SUBROUTINE MATINT_D(L,NR,ndim)
      USE NUMBERS
      USE PARINF_MODULE
      IMPLICIT REAL(FLOAT) (A-H,O-Z)
      DIMENSION L(Ndim,Ndim)
      DO M=1,NR,20
      K=MIN(M+19,NR)
      WRITE(IOUT,40)(J,J=M,K)
      DO I=1,NR
      WRITE(IOUT,50)I,(L(I,J),J=M,K)
      ENDDO
      ENDDO
      RETURN
40    FORMAT(/7X,20I6)
50    FORMAT(I4,3X,20(I6))
      END
      Program cryapi_inp
      Use numbers
      Use parame_module
      Use parinf_module
      Implicit real( float ) ( a-h, o-z )
      call cryread
      CALL KREDIN
      STOP
      END program cryapi_inp