 SUBROUTINE ReadInput_Vasp

! This routine reads an input template for Vasp

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

  CHARACTER(LCHARS) ::  Line
  INTEGER(KINT) ::  Iat
  INTEGER(KINT) :: I,J

  INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na
  INTEGER(KINT) :: NbType, NbTypeUser

  REAL(KREAL) :: B(3),xtmp,ytmp,ztmp

  LOGICAL :: TChk

  LOGICAL :: Debug


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

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

  FPBC=.TRUE.

 if (Input/="VASP") THEN

     ! Input was not Vasp, so many parameters are missing like lattice 
     ! constants...
     ! we read them now !
     ALLOCATE(FFF(3,nat))
     ! First geometry is a bit special for VASP as we have to set
     ! many things
     IF (DEBUG) WRITE(*,*) "Reading Vasp Parameters"
     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=998) NbAtType
998  CONTINUE
     NbType=0
     DO WHILE (NbAtType(NbType+1).NE.0)
        NbType=NbType+1
     END DO

     ! Do we know the atom types ?
     IF (AtTypes(1).EQ.'  ') THEN
        ! user has not provided atom types... we have to find them 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:2).NE.'US')
              READ(IOTMP,'(A)') Line
              Line=AdjustL(Line)
           END DO
           Line=adjustl(Line(3:))
           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)

! PFL 2011 Mar 8 ->
! We have to read the FFF flags :
     DO I=1,Nat
        READ(IOIN,*) Xtmp,ytmp,ztmp,FFF(1:3,I)
        DO J=1,3
           FFF(J,I)=AdjustL(FFF(J,I))
           CALL Upcase(FFF(J,I))
        END DO
     END DO
! <- PFL 2011 Mar 8 

  END IF

  ! In the case of VASP there is always the problem  of moving from one side
  ! of the box to the other...

     Renum=.TRUE.

     ! V_direct has been set in Read_geom
     IF (V_direct(1:6).EQ.'DIRECT') THEN
        Latr(1:3,1)=Lat_a
        Latr(1:3,2)=Lat_b
        Latr(1:3,3)=Lat_c
        B=1.
! TO DO: replace by LAPACK
        CALL Gaussj(Latr,3,3,B,1,1)
     ELSE
        Latr=0.
        Latr(1,1)=1.d0
        Latr(2,2)=1.d0
        Latr(3,3)=1.d0
     END IF

     ! Actualization of Frozen using the FFFF... 
     ! Frozen has the advantage ie if given, it imposes _ALL_ the FFF flags.
     IF (Frozen(1).NE.0) THEN
        WRITE(IOOUT,*) "Frozen is given. Flags of the given POSCAR are overriden"
        FFF='T'

        NFroz=0
        DO WHILE (Frozen(NFroz+1).NE.0)
           NFroz=NFroz+1
           FFF(1:3,Frozen(NFroz))='F'
        END DO
     ELSE
        WRITE(IOOUT,*) "Frozen not given : using  Flags of the given POSCAR"
        NFroz=0
        Frozen=0
        DO I=1,Nat
           IF ((FFF(1,I).EQ.'F').OR.(FFF(2,I).EQ.'F').OR.(FFF(3,I).EQ.'F')) THEN
              FFF(1:3,I)='F'
              NFroz=NFroz+1
              Frozen(NFroz)=I
           END IF
        END DO
        WRITE(IOOUT,*) "Frozen atoms:",Frozen(1:NFroz)
     END IF

     IF (Vmd) THEN
        if (debug) WRITE(*,*) "DBG MAIN, L803, VMD=T,NbTypes,AtTypes",NbTypes,AtTypes(1:NbTypes)
        Line=""
        DO I=1,NbTypes
           Line=TRIM(Line) // " " // TRIM(AdjustL(AtTypes(I))) 
        END DO
        Vasp_Title=Trim(Line) // " " // Trim(adjustL(Vasp_Title))
        if (debug) WRITE(*,*) "DBG MAIN, L809, VMD=T, Vasp_Title=",TRIM(Vasp_Title)
     END IF

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

   END SUBROUTINE READINPUT_VASP

