Chimie Théorique » Opt'n Path

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