SUBROUTINE ReadGeom_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 Path_module
  Use Io_module

  IMPLICIT NONE


  CHARACTER(LCHARS) :: LineTmp,Line

  INTEGER(KINT), ALLOCATABLE :: NbAtType(:) !na
  INTEGER(KINT) :: NbType, NbTypeUser
  INTEGER(KINT) :: I, J, Iat, NAtP, Idx
  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').or.valid('ReadGeom_vasp')

 if (debug) Call Header("Entering ReadGeom_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))
! PFL 2013.12.07 ->
! From VASP 5.3, POSCAR can have an additional line containing
! the names of the atoms
     VASP5=.FALSE.
     READ(IOIN,'(A)') Vasp_types
     LineTMp=AdjustL(Trim(Vasp_types))
     if (ICHAR(LineTmp(1:1))<48) THEN
        Call Warning("ReadGeom_Vasp","Illegal caracter found in reading Vasp types")
        STOP
     END IF
     If (ICHAR(LineTmp(1:1))>58) THEN
! we have a letter, the user has surely provided the names.
! We follow VASP and do not use them.
        VASP5=.TRUE.
        Vasp_types_user=Vasp_types
! we read the usual "number of atom by type" line
        READ(IOIN,'(A)') Vasp_types
     END IF
! <- PFL 2013.12.07     

     ! 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
! PFL 2013.12.07 ->
! A comment: in fact, this line might just not be 'Selective Dynamics'
     READ(IOIN,'(A)') Vasp_comment
     LineTMp=AdjustL(Trim(Vasp_comment))
     Call upcase(LineTMP)
     if (LineTMP(1:1)=="S") THEN
        Vasp_SelectD=.TRUE.
        READ(IOIN,'(A)') Vasp_direct
     ELSE
        Vasp_SelectD=.FALSE.
        Vasp_direct=Vasp_comment
     END IF
! <- PFL 2013.12.07

     V_direct=adjustl(Vasp_direct)
     Call upcase(V_direct)

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

     DO I=1,Nat
        IF (Vasp_SelectD) THEN
           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
        ELSE
           READ(IOIN,*) XyzGeomI(IGeom,1:3,I)
           FFFl(1:3,i)=.TRUE.
           FFF(1:3,i)="T"
        END IF
     END DO
! PFL 2013.12.07 ->
! Reading the VASP manual more carefully, it is stated:
! ! The seventh line (or eighth line if 'selective dynamics' is switched on)
! ! specifies whether the atomic positions are provided in cartesian coordinates
! ! or in direct coordinates (respectively fractional coordinates). As in the
! ! file KPOINTS only the first character on the line is significant and the
! ! only key characters recognized by VASP are 'C', 'c', 'K' or 'k' for
! ! switching to the cartesian mode.    
!
! We thus change V_direct to Direct if its first letter is not 'C' or 'K' 
! (we have it all in upper case)
     IF ((V_direct(1:1).NE.'C').AND.(V_direct(1:1).NE.'K')) THEN
        V_direct='DIRECT'
     END IF
! <- PFL 2013.12.07

     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.
! 
!! PFL 2013 Dec 07 ->
!!!!! For now we assume that all "POSCAR" used follow the same
!!!! pattern as the first one
!!! this has to be changed: reading the VASP header should be in
!!! a subroutine by itself and combined afterwards.
!!!! <- PFL 2013 Dec 07

! 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)
! PFL 2013 Dec 07 ->
! From VASP 5.3, POSCAR can have an additional line containing
! the names of the atoms
! Even thought VASP5 flag has been set on the first geometry,
! we are quite flexible and let the user mix old and new POSCAR:
        READ(IOIN,'(A)') LineTmp
        LineTMp=AdjustL(Trim(Vasp_types))
        if (ICHAR(LineTmp(1:1))<48) THEN
           Call Warning("ReadGeom_Vasp","Illegal caracter found in reading Vasp types")
           STOP
        END IF
        If (ICHAR(LineTmp(1:1))>58) THEN
! we have a letter, the user has surely provided the names.
! we read the usual "number of atom by type" line
           READ(IOIN,*) LineTmp
        !        if (debug) WRITE(*,*) "6 LineTmp:",Trim(LineTmp)
        END IF
! Do we have Selective Dynamics ?
        READ(IOIN,'(A)') Line
        LineTMp=AdjustL(Trim(Line))
        Call upcase(LineTMP)
        if (LineTMP(1:1)=="S") THEN
           if (.NOT.Vasp_SelectD) Vasp_comment=AdjustL(Line)
           Vasp_SelectD=.TRUE.
           READ(IOIN,'(A)') Vdir_loc
        ELSE
           Vasp_SelectD=.FALSE.
           Vdir_loc=Line
     END IF
! <- PFL 2013.12.07


! V_direct
        Vdir_loc=adjustl(Vdir_loc)
        Call upcase(Vdir_loc)
! PFL 2013.12.07 ->
     IF ((Vdir_loc(1:1).NE.'C').AND.(Vdir_loc(1:1).NE.'K')) THEN
        Vdir_loc='DIRECT'
     END IF
! <- PFL 2013.12.07


        !        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
           IF (Vasp_SelectD) THEN
              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
           ELSE
              READ(IOIN,*) XyzGeomI(IGeom,1:3,I)
           END IF

        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

     DO I=1,Nat
        DO J=1,3
           if (FFFl(j,i)) FFF(j,i)="T"
        END DO
     END DO

 Deallocate(FFFl,FFFt)


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

END SUBROUTINE ReadGeom_vasp
