Chimie Théorique » Opt'n Path

 SUBROUTINE ReadInput_siesta

! This routine reads an input template for Siesta

!----------------------------------------------------------------------

!  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

    FUNCTION SearchInput(Input,String,Line,Clean)  Result (Found)

      Use Vartypes

      use io_module

! Input

      TYPE (Input_line), POINTER, INTENT(IN) :: Input

      CHARACTER(*), INTENT(IN) :: String

      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean

! Output

      TYPE (Input_line), POINTER, INTENT(OUT) :: Line

      LOGICAL :: Found

    END FUNCTION SearchInput

    FUNCTION InString(Line,String,Case,Clean,Back)  Result(Pos)

      Use VarTypes

      implicit none

! Input

      CHARACTER(*), INTENT(IN) :: Line

      CHARACTER(*), INTENT(IN) :: String

      LOGICAL, OPTIONAL, INTENT(IN) :: Case

      LOGICAL, OPTIONAL, INTENT(IN) :: Back

      CHARACTER(*), OPTIONAL, INTENT(IN) :: Clean

! Output

! the position of String in Line (the first one) unless Back is present

      INTEGER(KINT) :: Pos

    END FUNCTION InString

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

      Use VarTypes

      Use io_module

      implicit none

!--------------------------------------------------------------- Input Variables

      character(len=*), intent(in)           :: routine, msg

      character(len=*), intent(in), optional :: file

      integer(KINT), intent(in), optional      :: line, unit

    END SUBROUTINE die

    SUBROUTINE Warning(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 Warning

 SUBROUTINE WriteList(Input,Unit)

! This routine reads an input template for Siesta

  use VarTypes

  use Io_module

  IMPLICIT NONE

! Input

      TYPE (Input_line), POINTER, INTENT(IN) :: Input

      INTEGER(KINT), OPTIONAL, INTENT(IN) :: Unit

    END SUBROUTINE WriteList

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!

 END  INTERFACE

!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  CHARACTER(132) ::  Line,Line2

  INTEGER(KINT) ::  Idx,LineL

  INTEGER(KINT) ::  IAt,I

  INTEGER(KINT) :: IoRead, ITmp

  REAL(KREAL) :: B(3),Xtmp, Ytmp, Ztmp

  REAL(KREAL) :: Lata,Latb,Latc,Alpha,Beta,Gamma

  REAL(KREAL), ALLOCATABLE :: GeomCartTmp(:,:) ! (Nat,3)

  REAL(KREAL) :: LatLoc(3,3)

  LOGICAL :: Debug

  LOGICAL :: FSpecies, FCoord

  TYPE(Input_Line), POINTER :: Search,Bla

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

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

  ! We read the Siesta input file

     ALLOCATE(Siesta_Input)

     NULLIFY(Siesta_Input%next)

     NULLIFY(Siesta_Input%prev)

     Current => Siesta_Input

     FSpecies=.FALSE.

     FCoord=.FALSE.

     READ(IOIN,'(A)',iostat=IoRead) Line

     DO WHILE (IoRead==0)

        Line=AdjustL(Line)

        if (debug) WRITE(*,*) 'Line:', Line

        current%Line=TRIM(Line)

        ALLOCATE(current%next)

        Current%next%prev => Current

        Current => Current%next

        current%line="toto"

        Nullify(Current%next)

        READ(IOIN,'(A)',iostat=IoRead) Line

     END DO

! With this procedure, the last current is not valid and should be deleted.

       Bla => Current%prev

       nullify(Bla%next)

       deallocate(Current)

       current => Bla

     if (debug) THEN

        WRITE(*,*) 'Input just read'

        Call WriteList(Siesta_input,Unit=IOOUT)

     END IF

! We analyse the input

! We look for the NumberofAtoms

     If (SearchInput(Siesta_Input,"NUMBEROFATOMS",Search,Clean=".-_")) THEN

        Line=AdjustL(Search%line)

        Idx=Index(Line," ")

        Line2=Trim(AdjustL(Line(Idx+1:)))

        Read(Line2,*) IAt

        if (Iat/=Nat) THEN

           Call Die('ReadInput_siesta','Nat in FDF sample different from  Nat Path input', Unit=IOOUT)

        END IF

     ELSE

! There is no atom number defined !!!

           Call Warning('ReadInput_siesta','No NumberOfAtoms in FDF sample', Unit=IOOUT)

           WRITE(current%Line,'(1X,"NumberOfAtoms ",I5)') Nat

           ALLOCATE(current%next)

           Current%next%prev => Current

           Current => Current%next

           Nullify(Current%next)

        END IF

! We look for the NumberofSpecies

        IF (SearchInput(Siesta_Input,"NUMBEROFSPECIES",Search,Clean=".-_")) THEN

           Line=AdjustL(Search%line)

           Idx=Index(Line," ")

           Line2=Trim(AdjustL(Line(Idx+1:)))

           Read(Line2,*) Siesta_NbSpecies

           ALLOCATE(Siesta_SpeciesMass(Siesta_NbSpecies))

           ALLOCATE(Siesta_SpeciesName(Siesta_NbSpecies))

        END IF

! We look for SystemLabel

        If (SearchInput(Siesta_Input,"SYSTEMLABEL",Search,Clean=".-_")) THEN

           Line=AdjustL(Search%line)

           Idx=Index(Line," ")

           Siesta_Label=Trim(adjustl(Line(Idx+1:)))

        ELSE

           Siesta_label='siesta'

        END IF

! We look for the ChemicalSpeciesLabel block

        IF (SearchInput(Siesta_Input,"CHEMICALSPECIESLABEL",Search,Clean=".-_")) THEN

           if (.NOT.ALLOCATED(IdxSpecies)) ALLOCATE(IdxSpecies(NAt))

           I=0

           Search => Search%next

           DO WHILE (InString(Search%line,'ENDBLOCK',Case=.FALSE.,Clean=".-_")==0)

              I=I+1

              if (I>Siesta_NbSpecies) THEN

                 Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel','Found more line in this block than NbSpecies !')

              END IF

              Line=AdjustL(Search%Line)

              Read(Line,*) ITmp

              Idx=Index(Line,' ')

              Line=AdjustL(Line(Idx+1:))

              Read(Line,*) Idx

              Siesta_SpeciesMass(ITmp)=Idx

              Idx=Index(Line,' ')

              Line=AdjustL(Line(Idx+1:))

              Idx=Index(Line,' ')

              Siesta_SpeciesName(ITmp)=Line(1:Idx-1)

              Search => Search%next

           END DO

           IF (I/=Siesta_NbSpecies) Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel', &

                'Number of lines in this block different from NbSpecies')

        ELSE

           Call Die('ReadInput_siesta:Reading ChemicalSpeciesLabel', &

                'Block ChemicalSpeciesLabel is mandatory !')

        END IF

! We look for the AtomicCoordinatesAndAtomicSpecies  block

        IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESANDATOMICSPECIES",Search,Clean=".-_")) THEN

           ALLOCATE(Siesta_Paste(Nat))

           Current=>Search

           DO I=1,NAt

              Search => Search%next

              Read(Search%line,*) XTmp,YTmp,ZTmp,IdxSpecies(I)

! We give a name to this atom

              IF (AtName(I)/="") THEN

!                    Write(Line,'(A,I5,1X,A2,1X,A32)') 'AtName & SpeciesName for atom ',I,AtName(I), &

!                         TRIM(Siesta_SpeciesName(IdxSpecies(I)))

!                    WRITE(*,'(1X,A)') Line

                 If (InString(Atname(I),TRIM(Siesta_SpeciesName(IdxSpecies(I))))==0) THEN

                    Write(Line,'(A,I5,1X,A2,1X,A32)') 'AtName /= SpeciesName for atom ',I,AtName(I), &

                         TRIM(Siesta_SpeciesName(IdxSpecies(I)))

                    Call Die('Readinput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)

                 END IF

! We look for the Mass number also

                 IF (Atome(I)/=Siesta_SpeciesMass(IdxSpecies(I))) THEN

                    Write(Line,'(A,I5,1X,I5,1X,I5)') 'AtMass /= SpeciesMass for atom ',I,Atome(I), &

                         Siesta_SpeciesMass(IdxSpecies(I))

                    Call Die('ReadInput_siesta:Reading AtomicCoordinatesAndAtomicSpecies',Line,Unit=IOOUT)

                 END IF

              ELSE

                 AtName(I)=AdjustL(Trim(Siesta_SpeciesName(IdxSpecies(I))))

                 Atome(I)=Siesta_SpeciesMass(IdxSpecies(I))

              END IF

! we look for something else at the end of the line

! We save everything but the x,y,z, and species description              

              Line=AdjustL(search%line)

! We skip x

              Idx=Index(Line,' ')

              Line=AdjustL(Line(Idx+1:))

! we skip y

              Idx=Index(Line,' ')

              Line=AdjustL(Line(Idx+1:))

! we skip z

              Idx=Index(Line,' ')

              Line=AdjustL(Line(Idx+1:))

! we skip species

              Idx=Index(Line,' ')

              Line=AdjustL(Line(Idx+1:))

              Siesta_Paste(I)=TRIM(Line)

           END DO

           if (debug) THEN

              WRITE(*,*) 'Input before deletion of coord block'

              Call WriteList(Siesta_input)

           END IF

! We will now delete this block from our input sample as it will then be

! written directly by Opt'n Path

! Search%next point on %endblock

           search => search%next

           IF (ASSOCIATED(Search%next)) THEN

! we are not at the end of the input file

              if (ASSOCIATED(Current%Prev)) THEN

                 Search%next%prev => current%prev

                 current%prev%next => search%next

              ELSE

! the coordinate block is at the begining of the input

                 siesta_input => Search%next

                 Nullify(Siesta_Input%Prev)

              END IF

           ELSE

! the coordinate block is the last part of the input

              nullify(current%prev%next)

           END IF

           if (debug) THEN

              WRITE(*,*) 'Input after deletion of coord block 1'

              Call WriteList(Siesta_input)

           END IF

           DO I=1,Nat+2

              bla=>current

              current=> current%next

              deallocate(bla)

           END DO

           if (debug) THEN

              WRITE(*,*) 'Input after deletion of coord block 2'

              Call WriteList(Siesta_input)

           END IF

        ELSE

           IF (SearchInput(Siesta_Input,"ZMATRIX",Search,Clean=".-_")) THEN

              call Die('ReadInput_Siesta','For now, I need the full block'//&

                   'AtomicCoordinatesAndAtomicSpecies, ZMatrix not yet handled.')

           ELSE

              call Die('ReadInput_Siesta','For now, I need the full block' //&

              'AtomicCoordinatesAndAtomicSpecies.')

           END IF

        END IF

        IF (SearchInput(Siesta_Input,"LATTICE",Search,Clean=".-_")) THEN

! This is a pbc calculation

           FPBC=.TRUE.

           IPer=3

           Kabeg=-1

           kaEnd=1

           kbBeg=-1

           kbEnd=1

           kcBeg=-1

           kcEnd=1

         IF (SearchInput(Siesta_Input,"LATTICECONSTANT",Search,Clean=".-_")) THEN

           Line=Adjustl(Search%Line)

! We discar the label

           Idx=Index(Line,' ')

           Line=AdjustL(Line(Idx+1:))

! we read the value

           Read(Line,*) Siesta_LatticeConstant

! We discard the value

           Idx=Index(Line,' ')

           Line=AdjustL(Line(Idx+1:))

           LineL=Len_Trim(Line)

! We read the unit

           If (LineL>0) THEN

              Call UpCase(Line)

              SELECT CASE (Line) 

                CASE ('ANG')

                   Siesta_Lat_Unit=1.d0

                CASE ('BOHR')

                   Siesta_Lat_Unit=a0

              END SELECT

           ELSE

              Siesta_Lat_Unit=1.d0

           END IF

         ELSE

! We have no LatticeConstant but we need one for our calculation

            Siesta_LatticeConstant=1.

            Siesta_Lat_Unit=1.d0

         END IF

         IF (SearchInput(Siesta_Input,"LATTICEVECTORS",Search,Clean=".-_")) THEN

           search => search%next

           Line=Adjustl(Search%Line)

           Read(Line,*) Lat_a(1:3)

           search => search%next

           Line=Adjustl(Search%Line)

           Read(Line,*) Lat_b(1:3)

           search => search%next

           Line=Adjustl(Search%Line)

           Read(Line,*) Lat_c(1:3)

         ELSEIF (SearchInput(Siesta_Input,"LATTICEPARAMETERS",Search,Clean=".-_")) THEN

           search => search%next

           Line=Adjustl(Search%Line)

           Read(Line,*) Lata,Latb,Latc,Alpha,Beta,Gamma

            Lat_a=0.

            Lat_b=0.

            Lat_c=0.

            Lat_a(1)=Lata

            Lat_b(1)=LatB*cos(Gamma*Pi/180d0)

            Lat_b(2)=LatB*sin(Gamma*Pi/180d0)

            Lat_c(1)=Latc*cos(Beta*Pi/180.)

            Lat_c(2)=(Latb*Latc*cos(Alpha*Pi/180.)-Lat_b(1)*Lat_c(1))/Lat_b(2)

            Lat_c(3)=sqrt(Latc*Latc-Lat_c(1)**2-Lat_c(2)**2)

         ELSE

! We have a lattice constant but nothing else. We issue a warning just in case,

! and we take the defaut values 1. 1. 1. 90. 90. 90.

            Line="LatticeConstant found, but not LatticeVectors nor LatticeParameters. Taking 1. 1. 1. 90. 90. 90."

            Call Warning("ReadInput_siesta",Line)

            Lat_a=(/1.,0.,0./)

            Lat_b=(/0.,1.,0./)

            Lat_c=(/0.,0.,1./)

         END IF

         IF (debug) THEN 

            WRITE(*,*) "Lattice vectors before scaling with LatticeConstant"

            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_a",Lat_a

            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_b",Lat_b

            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_c",Lat_c

         END IF

           Lat_a=Lat_a*Siesta_LatticeConstant*Siesta_Lat_Unit

           Lat_b=Lat_b*Siesta_LatticeConstant*Siesta_Lat_Unit

           Lat_c=Lat_c*Siesta_LatticeConstant*Siesta_Lat_Unit

         IF (debug) THEN 

            WRITE(*,*) "Lattice vectors After scaling with LatticeConstant"

            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_a",Lat_a

            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_b",Lat_b

            WRITE(*,'(1X,A,3(1X,F10.6))')   "Lat_c",Lat_c

         END IF

      END IF

         IF (SearchInput(Siesta_Input,"ATOMICCOORDINATESFORMAT",Search,Clean=".-_")) THEN

           Line=Adjustl(Search%Line)

           Idx=Index(Line,' ')

           Line=AdjustL(Line(Idx+1:))

           Call UpCase(Line)

           SELECT CASE (Line) 

              CASE ('ANG','NOTSCALEDCARTESIANANG')

                 Siesta_Unit_Read=1.d0

              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')

                 IF (FPBC) THEN

                    Siesta_Unit_Read=1.d0

                    V_DIRECT="DIRECT"

                    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)

                    LatLoc(1:3,1)=Lat_a

                    LatLoc(1:3,2)=Lat_b

                    LatLoc(1:3,3)=Lat_c

                    Allocate(GeomCartTmp(Nat,3))

                    Do I=1,NGeomI

                       GeomCartTmp=Reshape(XyzGeomI(I,1:3,1:Nat),(/Nat,3/),ORDER=(/2,1/))

! TO DO: shall we replace MatMull by DGEMM ?

                       XyzGeomI(I,1:3,1:Nat)=Reshape(MatMul(GeomCartTmp,Latloc),(/3,Nat/),ORDER=(/2,1/))

                   END DO

                 ELSE

                    WRITE(Line2,*) "AtomicCoordinatesFormat=" // Trim(Line) &

               //", but there is no information about Lattice parameters"

                    Call Die("ReadInput_siesta",Line2)

                 END IF

              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')

                 Siesta_Unit_Read=a0

                 IF (INPUT=='SIESTA') THEN

! We have read the coordinates, but not the unit. This is corrected here

                    XyZGeomI=XyzGeomI*a0

                 END IF

              CASE ('SCALEDCARTESIAN')

                 Siesta_Unit_Read=Siesta_LatticeConstant

                 If (INPUT=='SIESTA') THEN

                    XyzGeomI=XyzGeomI*Siesta_Unit_Read

                 END IF

            END SELECT

         ELSE

            Siesta_Unit_Read=1.d0

         END IF

        IF (SearchInput(Siesta_Input,"ATOMICCOORFORMATOUT",Search,Clean=".-_")) THEN

           Line=Adjustl(Search%Line)

           Idx=Index(Line,' ')

           Line=AdjustL(Line(Idx+1:))

           Call UpCase(Line)

           SELECT CASE (Line) 

              CASE ('ANG','NOTSCALEDCARTESIANANG')

                 Siesta_Unit_Write=1.d0

              CASE ('FRACTIONAL','SCALEDBYLATTICEVECTORS')

                 Siesta_Unit_Write=1.d0

                 V_DIRECT_Write='DIRECT'

              CASE ('BOHR','NOTSCALEDCARTESIANBOHR')

                 Siesta_Unit_Write=Unit ! 1/a0

            END SELECT

         END IF

         IF (SearchInput(Siesta_Input,"SUPERCELL",Search,Clean=".-_")) THEN

            Line=Adjustl(search%line)

            Call Upcase(Line)

            Call CleanString(Line,".-_")

            If (Index(Line,"%BLOCK")/=0) THEN

               Call Die('ReadInput_siesta:SuperCell',"SuperCell  NOT possible in Opt'n Path")

            END IF

         END IF

!         Call Die('Readinput_siesta:fin',' Lecture finie')

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

   END SUBROUTINE READINPUT_SIESTA