 SUBROUTINE ReadAnaList
! This routines read a list of geometrical variables to monitor
! This is inspired from Xyz2Path (that was inspired by Xyz2scan ...)

!----------------------------------------------------------------------
!  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

  LOGICAL :: Debug
  INTEGER(KINT) :: At1,At2,At3,At4
  INTEGER(KINT) :: I,Idx,J

  CHARACTER(LCHARS) :: Line

  Debug=Valid('ReadAnaList')
  If (Debug) Call Header("Entering ReadAnaList")

  Allocate(Bary)
  Nullify(Bary%Next)
  CurBary => Bary

  ALLOCATE(GeomList)
  Nullify(GeomList%Next)
  CurVar => GeomList

  NbCom=0
  NbVar=0
  FormAna='(1X,F8.3'
  DO I=1, Nb
     READ(IOIN,'(A)') Line
     Line=AdjustL(Line)
     Call Upcase(Line)

     SELECT CASE (Line(1:1))
        CASE ('B') 
! this is a bond
           CurVar%Type="BOND"
           Idx=Index(Line," ")
           Line=Line(Idx+1:)
           Read(Line,*) At1, At2
           CurVar%At1=At1
           CurVar%At2=At2
           CurVar%Value=0.
!           CurVar%PrintFactor=1.
           Curvar%SignDihedral=1
           if (debug) THEN
              WRITE(*,'("# b ",2I3)') At1,At2
           END IF
           FormAna=TRIM(FormAna) //',1X,F7.3'
           Allocate(CurVar%Next)
           CurVar => CurVar%Next
           Nullify(CurVar%Next)
        CASE ('A') 
! this is a valence angle
           CurVar%Type="ANGLE"
           Idx=Index(Line," ")
           Line=Line(Idx+1:)
           Read(Line,*) At1, At2,At3
           CurVar%At1=At1
           CurVar%At2=At2
           CurVar%At3=At3
           CurVar%Value=0.
!           CurVar%PrintFactor=180./Pi
           Curvar%SignDihedral=1
           if (debug) THEN
              WRITE(*,'("# a ",4(I3))') At1,At2,At3
           END IF
           FormAna=TRIM(FormAna) //',1X,F7.2'
           Allocate(CurVar%Next)
           CurVar => CurVar%Next
           Nullify(CurVar%Next)
        CASE ('D') 
! this is a dihedral
           CurVar%Type="DIHEDRAL"
           Idx=Index(Line," ")
           Line=Line(Idx+1:)
           Read(Line,*) At1, At2,At3,At4
           CurVar%At1=At1
           CurVar%At2=At2
           CurVar%At3=At3
           CurVar%At4=At4
           CurVar%Value=0.
!           CurVar%PrintFactor=180./Pi
           Curvar%SignDihedral=1
           if (debug) THEN
              WRITE(*,'("# d ",4(I3))') At1,At2,At3,At4
           END IF
           FormAna=TRIM(FormAna) //',1X,F7.2'
           Allocate(CurVar%Next)
           CurVar => CurVar%Next
           Nullify(CurVar%Next)
        CASE ('C') 
           NbCom=NbCom+1
           Idx=Index(Line," ")
           Line=Line(Idx+1:)
           READ(Line,*) At1
           Allocate(CurBary%ListAtoms(0:At1))
           Allocate(CurBary%Weights(1:At1))
           CurBary%Weights=1.
           CurBary%ListAtoms(0)=At1
           Idx=Index(Line," ")
           Line=Line(Idx+1:)
           READ(Line,*) CurBary%ListAtoms(1:At1)
           if (debug) THEN
              WRITE(*,'("# c ",I4,20(I3))') At1, &
                   (CurBary%ListAtoms(j),j=1,At1)
           END IF
           Allocate(CurBary%Next)
           CurBary => CurBary%Next
           Nullify(CurBary%Next)

        CASE Default
           Call Die('ReadAnaList','Variable not recognized: ' //Line, Unit=IOOUT)
     END SELECT

  END DO

  NbVar=Nb-NbCom

  CurVar => GeomList
  I=0

  Allocate(PrintGeomFactor(NbVar))
  DO WHILE (associated(CurVar%Next))
     I=I+1
     SELECT CASE (CurVar%TYPE)
        CASE ('ANGLE','DIHEDRAL')
           PrintGeomFactor(I)=180./pi
        CASE DEFAULT
           PrintGeomFactor(I)=1.
     END SELECT
     CurVar => CurVar%Next
  END DO
  
  FormAna=Trim(FormAna) // ',1X,F7.2,1X,F15.6)'


!  if (debug) WRITE(*,*) "NbVar, NbCom,Nb",NbVar,NbCom,Nb

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

END SUBROUTINE ReadAnaList
