 SUBROUTINE ReadInput_gaussian

! This routine reads an input template for Gaussian

!----------------------------------------------------------------------
!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon, 
!  Centre National de la Recherche Scientifique,
!  Université Claude Bernard Lyon 1. All rights reserved.
!
!  This work is registered with the Agency for the Protection of Programs 
!  as IDDN.FR.001.100009.000.S.P.2014.000.30625
!
!  Authors: P. Fleurat-Lessard, P. Dayal
!  Contact: optnpath@gmail.com
!
! This file is part of "Opt'n Path".
!
!  "Opt'n Path" is free software: you can redistribute it and/or modify
!  it under the terms of the GNU Affero General Public License as
!  published by the Free Software Foundation, either version 3 of the License,
!  or (at your option) any later version.
!
!  "Opt'n Path" is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU Affero General Public License for more details.
!
!  You should have received a copy of the GNU Affero General Public License
!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.
!
! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,
! for commercial licensing opportunities.
!----------------------------------------------------------------------

  use VarTypes
  use Path_module
  use Io_module

  IMPLICIT NONE

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


    SUBROUTINE die(routine, msg, file, line, unit)

      Use VarTypes
      Use io_module

      implicit none
      character(len=*), intent(in)           :: routine, msg
      character(len=*), intent(in), optional :: file
      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE die


  END INTERFACE


  CHARACTER(132) ::  Line,LineUp
  INTEGER(KINT) :: LineL, Idx, Iat
  INTEGER(KINT) :: I,NTmp
  REAL(KREAL) :: Lat(3,3)
  
  LOGICAL :: Debug


  Debug=Valid("readinput").OR.Valid("readinput_gaussian")

 if (debug) Call Header("Entering ReadInput_Gaussian")

  ! We read the Gaussian input file
  ! First, the root
     IF (DEBUG) WRITE(*,*) "Reading Gauss Root"
     ALLOCATE(Gauss_Root)
     NULLIFY(Gauss_Root%next)
     Current => Gauss_root
     LineL=1
     DO WHILE (LineL.NE.0)
        READ(IOIN,'(A)') Line
        Line=AdjustL(Line)
        LineL=len_Trim(Line)
! we do not read the %chk line
        Idx=INDEX(Line,"chk")
        IF ((LineL.NE.0).AND.(Idx.EQ.0)) THEN
           current%Line=TRIM(Line)
           ALLOCATE(current%next)
           Current => Current%next
           Nullify(Current%next)
        END IF
     END DO

!     Current => Gauss_root
!     DO WHILE (ASSOCIATED(Current%next))
!        WRITE(*,'(1X,A)') Trim(current%line)
!        Current => current%next
!     END DO

     ! Now the comment... 
     IF (DEBUG) WRITE(*,*) "Reading Gauss Comment"
     ALLOCATE(Gauss_Comment)
     NuLLIFY(Gauss_Comment%Next)
     Current => Gauss_comment
     LineL=1
     DO WHILE (LineL.NE.0)
        READ(IOIN,'(A)') Line
        Line=AdjustL(Line)
        LineL=len(Trim(Line))
        IF (LineL.NE.0) THEN
           current%Line=TRIM(Line)
           ALLOCATE(current%next)
           Current => Current%next
           Nullify(Current%next)
        END IF
     END DO

 !    Current => Gauss_comment
 !    DO WHILE (ASSOCIATED(Current%next))
 !       WRITE(*,'(1X,A)') Trim(current%line)
 !       Current => current%next
 !    END DO

     ! Now the charge
     IF (DEBUG) WRITE(*,*) "Reading Gauss Charge"
     READ(IOIN,'(A)') Gauss_Charge
     if (debug) WRITE(*,*) "Gauss_charge=",TRIM(Gauss_charge)
     ! We now read the Paste part...
     ALLOCATE(Gauss_Paste(NAt))
     LineL=1
     Iat=0
     IPer=0
     FPBC=.FALSE.
     Gauss_paste=" "
     DO While (LineL.GT.0)
        READ(IOIN,'(A)') Line
        Line=AdjustL(Line)
        LineL=Len_TRIM(Line)
        IF (LineL.GT.0) THEN
           LineUp=Line
           Call UpCase(LineUp)
           if (LineUp(1:2)=="TV") THEN
              FPBC=.TRUE.
              IPer=IPer+1
              If (Iper>3) THEN
                 Call Die("ReadInput Gaussian","Iper>3",Unit=IOOUT)
              END IF
              NTmp=Index(LineUp," ")
              LineUp=LineUp(NTmp:)
              Read(LineUp,*) Lat(IPer,1:3)
           ELSE
              Iat=Iat+1
! we search for additional information at the end of the line
! for example ONIOM layers
! TO detect the end of the line, we use the fact that all reals
! should contain a '.', and that we have 3 reals/line.
              Idx=Index(Line,'.')
              Line=ADJUSTL(Line(Idx+1:))
              Idx=Index(Line,'.')
              Line=ADJUSTL(Line(Idx+1:))
              Idx=Index(Line,'.')
              Line=ADJUSTL(Line(Idx+1:))
              Idx=Index(Line,' ')
              If (Idx>0) THEN
                 Line=ADJUSTL(Line(Idx:))
                 IF (LEN_TRIM(Line).GT.0) THEN 
                    Gauss_paste(Iat)=ADJUSTL(TRIM(Line))
                 END IF
              ELSE
                 Gauss_paste(Iat)=""
              END IF
           END IF
        END IF
           
     END DO


     IF (Iat.NE.Nat) THEN
        WRITE(Line,*) "I found ", Iat," lines for the geometry instead of ",Nat
        Call Die("ReadInput Gaussian","Line",UNIT=IOOUT)
     END IF

  IF (FPBC) THEN
     Lat_a(1:3)=Lat(1,1:3)
     Lat_b(1:3)=Lat(2,1:3)
     Lat_c(1:3)=Lat(3,1:3)
     If (IPer>=1) THEN
        kaBeg=-1
        kaEnd=1
     END IF
     If (IPer>=2) THEN
        kbBeg=-1
        kbEnd=1
     END IF
     If (IPer==3) THEN
        kcBeg=-1
        kcEnd=1
     END IF
     If (IPer>3) THEN
        Call Die("Readinput_gaussian","Found too many Tv lines !",Unit=IOOUT)
     END IF
  END IF


     ! We now read the last part
     IF (DEBUG) WRITE(*,*) "Reading Gauss End"
     !     READ(IOIN,'(A)') Line
     ALLOCATE(Gauss_End)
     NuLLIFY(Gauss_End%Next)
     Current => Gauss_End
     LineL=1
     DO WHILE (1.EQ.1)
        READ(IOIN,'(A)',END=999) Line
        Line=AdjustL(Line)
        LineL=len(Trim(Line))
        current%Line=TRIM(Line)
        ALLOCATE(current%next)
        Current => Current%next
        Nullify(Current%next)
     END DO
999  CONTINUE

     IF (Debug) THEN
        ! Write the gaussian input file for testing purposes
        Current => Gauss_root
        DO WHILE (ASSOCIATED(Current%next))
           WRITE(*,'(1X,A)') Trim(current%line)
           Current => current%next
        END DO
        
        WRITE(*,*) 
!        WRITE(*,*) '//INFO// Comment original:'
        
        Current => Gauss_comment
        DO WHILE (ASSOCIATED(Current%next))
           WRITE(*,'(1X,A)') Trim(current%line)
           Current => current%next
        END DO
        
        WRITE(*,*) 
        WRITE(*,*) Trim(Gauss_charge)
        
        DO I=1,Nat
           WRITE(*,'(1X,A10,3(1X,F15.8),1X,A)') Trim(AtName(I)),XyzGeomI(1,1:3,I), TRIM(Gauss_Paste(I))
        END DO
        
        WRITE(*,*) 
        Current => Gauss_End
        DO WHILE (ASSOCIATED(Current%next))
           WRITE(*,'(1X,A)') Trim(current%line)
           Current => current%next
        END DO
        
        WRITE(*,*) 

        Call Header("Exiting ReadInput_Gaussian")

     END IF



   END SUBROUTINE READINPUT_GAUSSIAN
