 SUBROUTINE ReadInput_Mopac

! This routine reads an input template for MOPAC

!----------------------------------------------------------------------
!  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(LCHARS) ::  Line,LineUp
  INTEGER(KINT) :: LineL, Idx, NTmp
  INTEGER(KINT) :: NatMopac
  REAL(KREAL) :: Lat(3,3)
  
  LOGICAL :: Debug


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

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

! The structure is:
! A MOPAC data set normally consists of one line of keywords, two lines of user-defined text, then the coordinates
! Then  a blank line or a line of 0.
! then the symmetry description.
! comment lines start with * and can be anywhere !!!

  ! First, the root
  IF (DEBUG) WRITE(*,*) "Reading Mopac input"
  ALLOCATE(Mopac_Root)
  NULLIFY(Mopac_Root%next)
  ALLOCATE(Mopac_Comment)
  NULLIFY(Mopac_Comment%next)
  ALLOCATE(Mopac_End)
  NuLLIFY(Mopac_End%Next)
  Current => Mopac_root
  CurCom => Mopac_Comment
  LineL=1
  NTmp=0
  DO WHILE (NTmp.LT.3)
     READ(IOIN,'(A)') Line
     Line=AdjustL(Line)
     LineL=len(Trim(Line))
     IF (Line(1:1)/="*") THEN
        IF (NTmp==0) THEN
           LineUp=Line
           Call UpCase(LineUp)
           Idx=Index(LineUp,'GRADIENTS')
           If (Idx==0) Line=TRIM(Line) // " GRADIENTS"
           Idx=Index(LineUp,'1SCF')
           If (Idx==0) Line=TRIM(Line) // " 1SCF"
        END IF
        current%Line=TRIM(Line)
        ALLOCATE(current%next)
        Current => Current%next
        Nullify(Current%next)
        NTmp=NTmp+1
     ELSE
        CurCom%Line=TRIM(LINE)
        ALLOCATE(CurCom%Next)
        CurCom => CurCom%Next
        NULLIFY(CurCom%Next)
     END IF
  END DO

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

! Now the geometry... that we just skip :)
! PFL 2013 Apr 
! We take care that there is no Translation vectors...
! We also check that the number of atoms is ok
  IF (DEBUG) WRITE(*,*) "Reading Mopac Geometry"
  Mopac_EndGeom=""
  LineL=1
  NatMopac=0
  Lat=0.d0
  IPer=0
  FPBC=.FALSE.
  DO WHILE (LineL.NE.0)
     READ(IOIN,'(A)',END=989) Line
     Line=AdjustL(Line)
     LineL=len(Trim(Line))
     ! The last line might be either blank or filled with 0
     If (LineL>0) THEN
        SELECT CASE (Line(1:1))
          CASE ("0") 
             LineL=0
             Mopac_EndGeom=Trim(Line)
          CASE("*") 
             CurCom%Line=TRIM(LINE)
             ALLOCATE(CurCom%Next)
             CurCom => CurCom%Next
             NULLIFY(CurCom%Next)
          CASE DEFAULT
             LineUp=Line
             Call UpCase(LineUp)
             If (LineUp(1:2)=="TV") THEN
                FPBC=.TRUE.
                IPer=IPer+1
                If (Iper>3) THEN
                   Call Die("ReadInput Mopac","Iper>3",Unit=IOOUT)
                END IF
                NTmp=Index(LineUp," ")
                LineUp=LineUp(NTmp:)
                Read(LineUp,*) Lat(IPer,1:3)
             ELSE
                NatMopac=NatMopac+1
             END IF
         END SELECT
      END IF
     
   END DO

!  WRITE(*,*) "NatMopac,Nat:",NAtMopac,Nat
  IF (NatMopac/=Nat) Call Die("ReadInput_mopac","Nat read does not mat nat",Unit=IOOUT)
  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_mopac","Found too many Tv lines !",Unit=IOOUT)
     END IF
  END IF

! If we are here, there might be something else to read: Mopac_end

  ! We now read the last part
  IF (DEBUG) WRITE(*,*) "Reading Mopac End"
  !     READ(IOIN,'(A)') Line
  Current => Mopac_End
  LineL=1
  DO WHILE (1.EQ.1)
     READ(IOIN,'(A)',END=989) Line
     Line=AdjustL(Line)
     LineL=len(Trim(Line))
     IF (Line(1:1)/="*") THEN
        current%Line=TRIM(Line)
        ALLOCATE(current%next)
        Current => Current%next
        Nullify(Current%next)
        NTmp=NTmp+1
     ELSE
        CurCom%Line=TRIM(LINE)
        ALLOCATE(CurCom%Next)
        CurCom => CurCom%Next
        NULLIFY(CurCom%Next)
        END IF
     END DO
989  CONTINUE

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

END SUBROUTINE READINPUT_Mopac
