Chimie Théorique » OpenPath

SUBROUTINE Read_Geom(input)

  Use Path_module

  Use Io_module

  IMPLICIT NONE

  CHARACTER(32), INTENT(IN) :: input    

  CHARACTER(132) :: LineTmp,Line

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

  INTEGER(KINT) :: NbType, NbTypeUser

  INTEGER(KINT) :: I, J, Iat,NAtP,Idx,JStart

  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')

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

  if (debug) WRITE(*,*) "Input:",Trim(Input)

  SELECT CASE(Input)

  CASE ('XYZ','CART')

     DO I=1,NGeomI

        IF (DEBUG) WRITE(*,*) "Reading Geom :",I

        READ(IOIN,*) NAtp

        if (NAtp.NE.Nat) THEN

           IF (I==1) 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))

           ELSE

              WRITE(IOOUT,*) 'Number of atoms not consistent between geometries. STOP'

              STOP

           END IF

        END IF

        READ(IOIN,'(A)') Line

        DO J=1,NAt

           READ(IOIN,*) AtName(J),XyzGeomI(I,1:3,J)

        END DO

        If (Debug) THEN

           WRITE(*,*) "Geom ",I

           DO J=1,NAt

              WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(J),XyzGeomI(I,1:3,J)

           END DO

        END IF

     END DO

  CASE ('TURBOMOLE')

     DO I=1,NGeomI

        IF (DEBUG) WRITE(*,*) "Reading Geom :",I

        JStart=1

! read the $coord line, or $user-defined, or $end

        Line='$ok'

        DO WHILE (Line(1:1).EQ."$")

           READ(IOIN,'(A)') Line 

           Line=AdjustL(Line)

        END DO

        J=1

        READ(Line,*) XyzGeomI(I,1:3,J),AtName(J)

        DO J=2,NAt

           READ(IOIN,*) XyzGeomI(I,1:3,J),AtName(J)

        END DO

        XyzGeomI(I,:,:)= XyzGeomI(I,:,:)*a0

        If (Debug) THEN

           WRITE(*,*) "Geom ",I

           DO J=1,NAt

              WRITE(*,'(1X,A2,3(1X,F15.6))') AtName(J),XyzGeomI(I,1:3,J)

           END DO

        END IF

     END DO

  CASE ('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))

     READ(IOIN,'(A)') Vasp_types

     ! 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

     READ(IOIN,'(A)') Vasp_comment

     READ(IOIN,'(A)') Vasp_direct

     V_direct=adjustl(Vasp_direct)

     Call upcase(V_direct)

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

     DO I=1,Nat

        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')

           END DO

     END DO

     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.

! 

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

! Attypes

        READ(IOIN,*) LineTmp

        !        if (debug) WRITE(*,*) "6 LineTmp:",Trim(LineTmp)

! Comment

        READ(IOIN,*) LineTmp,line

        !        if (debug) WRITE(*,*) "7 LineTmp:",Trim(LineTmp),Trim(Line)

! V_direct

        READ(IOIN,*) Vdir_loc

        Vdir_loc=adjustl(Vdir_loc)

        Call upcase(Vdir_loc)

        !        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

           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

        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

! PFL 23 Nov 2008 ->

! This is done in Path.f90

! ! We put FFrozen in place

!   Idx=1

!   DO I=1,Nat

!      TChk=.FALSE.

!      DO J=1,3

!         TChk=Tchk.OR.(.NOT.FFFl(j,i))

!      END DO

!      IF (TChk) THEN

!         Frozen(Idx)=I

!         Idx=Idx+1

!      END IF

!   END DO

!   Frozen(Idx)=0

! We put FFF in place because FFFl is a logical flag

     DO I=1,Nat

        DO J=1,3

           if (FFFl(j,i)) FFF(j,i)="T"

        END DO

     END DO

!<- PFL 23 Nov 2008

 Deallocate(FFFl,FFFt)

  CASE Default

     WRITe(*,*) 'Input=',trim(Input),' UNKNOWN. Stop'

     STOP

  END SELECT

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

END SUBROUTINE Read_Geom