SUBROUTINE Read_Geom(input)

  Use Path_module
  Use Io_module

  IMPLICIT NONE


  CHARACTER(32), INTENT(IN) :: input    
  CHARACTER(132) :: LineTmp,Line

  INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na
  INTEGER(KINT) :: NbType, NbTypeUser
  INTEGER(KINT) :: I, J, Iat,NAtP,Idx,JStart
  INTEGER(KINT) :: IGeom
  REAL(KREAL) :: Xt,Yt,Zt
  LOGICAL :: Debug, TChk
  CHARACTER(4), ALLOCATABLE :: FFFt(:,:) ! 3,Na
  LOGICAL, ALLOCATABLE :: FFFl(:,:) ! 3,Na

!!! for VASP inputs
! latx_loc used to check that the lattice parameters are the same for all geometries
  REAL(KREAL) :: lata_loc(3), latb_loc(3), latc_loc(3)
! VaspPar_loc is the local value of the vasp parametre
  REAL(KREAL) :: VaspPar_loc
! Vdir_loc is used to convert the read geometry to Cartesian.
! V_direct is set by the first geometry.
  CHARACTER(LCHARS) :: Vdir_loc


  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID
  END INTERFACE

  debug=valid('Read_geom')
 if (debug) Call Header("Entering Read_Geom")
  if (debug) WRITE(*,*) "Input:",Trim(Input)
  if (debug) WRITE(*,*) "NgeomI:",NGeomI

  SELECT CASE(Input)
  CASE ('XYZ','CART')
     DO I=1,NGeomI
        IF (DEBUG) WRITE(*,*) "Reading Geom :",I
        READ(IOIN,*) NAtp
        if (NAtp.NE.Nat) THEN
           IF (I==1) THEN
              WRITE(IOOUT,*) 'WARNING Number of atoms not consistent between NAMELIST &path and First geom'
              WRITE(IOOUT,*) "Using:",Natp 
              DEALLOCATE(XyzGeomI, AtName)
              Nat=Natp 
              ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
           ELSE
              WRITE(IOOUT,*) 'Number of atoms not consistent between geometries. STOP'
              STOP
           END IF
        END IF
        READ(IOIN,'(A)') Line
        DO J=1,NAt
           READ(IOIN,*) AtName(J),XyzGeomI(I,1:3,J)
        END DO
        If (Debug) THEN
           WRITE(*,*) "Geom ",I
           DO J=1,NAt
              WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(J),XyzGeomI(I,1:3,J)
           END DO
        END IF
     END DO
  CASE ('TURBOMOLE')
     DO I=1,NGeomI
        IF (DEBUG) WRITE(*,*) "Reading Geom :",I
        JStart=1
! read the $coord line, or $user-defined, or $end
        Line='$ok'
        DO WHILE (Line(1:1).EQ."$")
           READ(IOIN,'(A)') Line 
           Line=AdjustL(Line)
        END DO
        J=1
        READ(Line,*) XyzGeomI(I,1:3,J),AtName(J)
        DO J=2,NAt
           READ(IOIN,*) XyzGeomI(I,1:3,J),AtName(J)
        END DO
        XyzGeomI(I,:,:)= XyzGeomI(I,:,:)*a0
        If (Debug) THEN
           WRITE(*,*) "Geom ",I
           DO J=1,NAt
              WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(J),XyzGeomI(I,1:3,J)
           END DO
        END IF
     END DO
  CASE ('VASP')
     ALLOCATE(FFF(3,nat),FFFl(3,Nat),FFFt(3,Nat))
     ! First geometry is a bit special for VASP as we have to set
     ! many things
     IGeom=1
     IF (DEBUG) WRITE(*,*) "Reading Geom :",IGeom
     READ(IOIN,'(A)') Vasp_Title
     READ(IOIN,*) Vasp_param

     READ(IOIN,*) Lat_a
     READ(IOIN,*) Lat_b
     READ(IOIN,*) Lat_c

     Lat_a=Lat_a*Vasp_param
     Lat_b=Lat_b*Vasp_param
     Lat_c=Lat_c*Vasp_param

     ALLOCATE(NbAtType(nat))
     READ(IOIN,'(A)') Vasp_types
     ! First, how many different types ?
     NbAtType=0
     READ(Vasp_types,*,END=999) NbAtType
999  CONTINUE
     NbType=0
     NatP=0
     DO WHILE (NbAtType(NbType+1).NE.0)
        NbType=NbType+1
        Natp=Natp+NbAtType(NbType)
     END DO


     if (NAtp.NE.Nat) THEN
        WRITE(IOOUT,*) 'WARNING Number of atoms not consistent between NAMELIST &path and First geom'
        WRITE(IOOUT,*) "Using:",Natp 
        DEALLOCATE(XyzGeomI, AtName)
        Nat=Natp 
        ALLOCATE(XyzGeomI(NGeomI,3,Nat), AtName(NAt))
     END if

     ! Do we know the atom types ?
     IF (AtTypes(1).EQ.'  ') THEN
        ! user has not provided atom types... we have to find them by ourselves
        ! by looking into the POTCAR file...
        INQUIRE(File="POTCAR", EXIST=Tchk)
        IF (.NOT.Tchk) THEN
           WRITE(*,*) "ERROR! No AtTypes provided, and POTCAR file not found"
           STOP
        END IF
        OPEN(IOTMP,File="POTCAR")
        DO I=1,NbType
           Line='Empty'
           DO WHILE (Line(1:5).NE.'TITEL')
              READ(IOTMP,'(A)') Line
              Line=AdjustL(Line)
           END DO
           Line=adjustl(Line(9:))
           Idx=index(Line," ")
           Line=adjustl(Line(Idx+1:))
           AtTypes(I)=Line(1:2)
        END DO
        if (debug) WRITE(*,'(A,100(1X,A2))') "ReadG:VASP AtTypes",AtTypes(1:NbType)
        CLOSE(IOTMP)

     ELSE  !AtTypes(1).EQ.'  '
        ! user has provided atom types
        NbTypeUser=0
        DO WHILE (AtTypes(NbTypeUser+1).NE.'  ')
           NbTypeUser=NbTypeUser+1
        END DO
        IF (NbType.NE.NbTypeUser) THEN
           WRITE(*,*) "ERROR Read_Geom : NbType in POSCAR do not match AtTypes"
           STOP
        END IF
     END IF

     IAt=1
     DO I=1,NbType
        DO J=1,NbAtType(I)
           AtName(Iat)=AtTypes(I)
           Iat=Iat+1
        END DO
     END DO
     DEALLOCATE(NbAtType)

     NbTypes=NbType

     READ(IOIN,'(A)') Vasp_comment
     READ(IOIN,'(A)') Vasp_direct

     V_direct=adjustl(Vasp_direct)
     Call upcase(V_direct)

     if (debug) WRITE(*,*) "V_direct=",Trim(V_Direct)

     DO I=1,Nat
        READ(IOIN,*) XyzGeomI(IGeom,1:3,I),FFF(1:3,I)
           DO J=1,3
              FFF(J,I)=AdjustL(FFF(J,I))
              FFFl(j,i)=(FFF(j,i)(1:1).EQ.'T').OR.(FFF(j,i)(1:1).EQ.'t')
           END DO

     END DO
     IF (V_direct(1:6).EQ.'DIRECT') THEN
        DO I=1,Nat
           Xt=XyzGeomI(IGeom,1,I)
           Yt=XyzGeomI(IGeom,2,I)
           Zt=XyzGeomI(IGeom,3,I)
           XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1)
           XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2)
           XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3)
        END DO

     END IF

     ALLOCATE(X0_vasp(Nat),Y0_vasp(Nat),Z0_vasp(Nat))
     X0_vasp=XyzGeomI(IGeom,1,:)
     Y0_vasp=XyzGeomI(IGeom,2,:)
     Z0_vasp=XyzGeomI(IGeom,3,:)

     if (debug)  THEN
        WRITE(*,*) Nat
        WRITE(*,*) "Geom ",IGeom,"/",NGeomI
        DO I=1,Nat
           WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(I),XyzGeomI(IGeom,1:3,I)
        END DO
     END IF


     DO IGeom=2,NGeomI
!!! PFL 16th March 2008 ->
! Users might want to use different conditions for different geometries
! i.e. to mix Cartesian and Direct :-p
! thus we check for each geom which condition is used.
! 

! Title
        READ(IOIN,*) LineTmp
        !        if (debug) WRITE(*,*) "1 LineTmp:",Trim(LineTmp)
! Vasp_Param
        READ(IOIN,*) VaspPar_loc
        !        if (debug) WRITE(*,*) "2 LineTmp:",Trim(LineTmp)
! Lattice parameters
        READ(IOIN,*) Lata_loc 
        !        if (debug) WRITE(*,*) "3 LineTmp:",Trim(LineTmp)
        READ(IOIN,*) latb_loc
        !        if (debug) WRITE(*,*) "4 LineTmp:",Trim(LineTmp)
        READ(IOIN,*) Latc_loc
        !        if (debug) WRITE(*,*) "5 LineTmp:",Trim(LineTmp)
! Attypes
        READ(IOIN,*) LineTmp
        !        if (debug) WRITE(*,*) "6 LineTmp:",Trim(LineTmp)
! Comment
        READ(IOIN,*) LineTmp,line
        !        if (debug) WRITE(*,*) "7 LineTmp:",Trim(LineTmp),Trim(Line)
! V_direct
        READ(IOIN,*) Vdir_loc
        Vdir_loc=adjustl(Vdir_loc)
        Call upcase(Vdir_loc)
        !        if (debug) WRITE(*,*) "8 LineTmp:",Trim(LineTmp)

        Lata_loc=Lata_loc*VaspPar_loc-Lat_a
        if (dot_product(lata_loc,Lata_loc).GE.1e-6) THEN
           WRITE(*,*) "Lattice parameter a is different between geometries 1 and ",IGeom
           WRITE(IOOUT,*) "Lattice parameter a is different between geometries 1 and ",IGeom
           STOP
        END IF

        Latb_loc=Latb_loc*VaspPar_loc-Lat_b
        if (dot_product(latb_loc,Latb_loc).GE.1e-6) THEN
           WRITE(*,*) "Lattice parameter b is different between geometries 1 and ",IGeom
           WRITE(IOOUT,*) "Lattice parameter b is different between geometries 1 and ",IGeom
           STOP
        END IF

        Latc_loc=Latc_loc*VaspPar_loc-Lat_c
        if (dot_product(latc_loc,Latc_loc).GE.1e-6) THEN
           WRITE(*,*) "Lattice parameter c is different between geometries 1 and ",IGeom
           WRITE(IOOUT,*) "Lattice parameter c is different between geometries 1 and ",IGeom
           STOP
        END IF

        DO I=1,Nat
           READ(IOIN,*) XyzGeomI(IGeom,1:3,I),FFFt(1:3,I)
           DO J=1,3
! A frozen atom in any geom is considered frozen in all geom
              FFFt(j,I)=AdjustL(FFFt(j,i))
              FFFl(j,i)=FFFl(j,i).AND.(FFFt(j,i).EQ.'T')
           END DO
        END DO

        IF (Vdir_loc(1:6).EQ.'DIRECT') THEN
           DO I=1,Nat
              Xt=XyzGeomI(IGeom,1,I)
              Yt=XyzGeomI(IGeom,2,I)
              Zt=XyzGeomI(IGeom,3,I)
              XyzGeomI(IGeom,1,I)=Xt*lat_a(1)+Yt*lat_b(1)+Zt*lat_c(1)
              XyzGeomI(IGeom,2,I)=Xt*lat_a(2)+Yt*lat_b(2)+Zt*lat_c(2)
              XyzGeomI(IGeom,3,I)=Xt*lat_a(3)+Yt*lat_b(3)+Zt*lat_c(3)
           END DO
        END IF

        if (debug)  THEN
           WRITE(*,*) Nat
           WRITE(*,*) "Geom ",IGeom,"/",NGeomI
           DO I=1,Nat
              WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(I),XyzGeomI(IGeom,1:3,I)
           END DO
        END IF

     END DO

! PFL 23 Nov 2008 ->
! This is done in Path.f90
! ! We put FFrozen in place
!   Idx=1
!   DO I=1,Nat
!      TChk=.FALSE.
!      DO J=1,3
!         TChk=Tchk.OR.(.NOT.FFFl(j,i))
!      END DO
!      IF (TChk) THEN
!         Frozen(Idx)=I
!         Idx=Idx+1
!      END IF
!   END DO
!   Frozen(Idx)=0
! We put FFF in place because FFFl is a logical flag
     DO I=1,Nat
        DO J=1,3
           if (FFFl(j,i)) FFF(j,i)="T"
        END DO
     END DO

!<- PFL 23 Nov 2008

 Deallocate(FFFl,FFFt)

  CASE Default
     WRITe(*,*) 'Input=',trim(Input),' UNKNOWN. Stop'
     STOP

  END SELECT


 if (debug) Call Header("Exiting Read_Geom")

END SUBROUTINE Read_Geom
