Chimie Théorique » OpenPath

SUBROUTINE Write_path(Iopt)

  Use Path_module

  Use Io_module

  IMPLICIT NONE

  INTEGER(KINT), INTENT(IN) :: Iopt

  INTEGER(KINT) :: IGeom, Idx, Iat, I,J

  REAL(KREAL), ALLOCATABLE :: GeomTmpC(:), GeomTmp(:),GeomTmpC2(:,:)

  CHARACTER(SCHARS) :: Line

  CHARACTER(LCHARS) ::Title

  LOGICAL :: Debug

  INTERFACE

     function valid(string) result (isValid)

       CHARACTER(*), intent(in) :: string

       logical                  :: isValid

     END function VALID

  END INTERFACE

  debug=valid('write_path').OR.valid('writepath')

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

  ALLOCATE(GeomTmpC(3*Nat), GeomTmp(NCoord),GeomTmpC2(Nat,3))

  IF (IOpt.GE.0) THEN

     WRITE(Line,'(I5)') Iopt

  ELSE

     Line="Ini"

     Energies=0.

  END IF

  OPEN(IOTMP,File=Trim(PathName) // '.' // AdjustL(TRIM(Line)))

  IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'MIXED').OR.(COORD.EQ.'BAKER')) &

       OPEN(IOCART,File=TRIM(PathName) // '_cart.' //  AdjustL(TRIM(Line)))

  DO IGeom=1,NGeomF

     WRITE(Title,"('Geometry ',I3,'/',I3,' for iteration ',I3,' E=',F13.6)") Igeom,NgeomF,Iopt,Energies(IGeom)

     SELECT CASE(Coord)

     CASE ('CART','HYBRID')

        ! XyzGeomF is 3,Nat... but Printgeom expects Nat,3... 

        ! this should really be rewritten !!!!

        GeomTmp=Reshape(reshape(XyzGeomF(IGeom,:,:),(/Nat,3/),ORDER=(/2,1/)),(/3*Nat/)) 

     CASE ('ZMAT','MIXED','BAKER')

        GeomTmp=IntCoordF(IGeom,:)

     CASE DEFAULT

        WRITE(*,*) "Coord=",TRIM(Coord)," not recognized in Write_path. L34.STOP"

        STOP

     END SELECT

     ! Nothing for the baker case.

     Call PrintGeom(Title,Nat,NCoord,GeomTmp,Coord,IoTmp,Atome,Order,OrderInv,IndZmat)

     IF (COORD.EQ.'ZMAT') THEN

        ! Writing the Cartesian coordinates

        WRITE(IOCART,'(1X,I5)') Nat

        WRITE(IOCART,'(1X,A)') Title

        ! we have to generate the cartesian coordinates from the internal coordinates

        Call Int2Cart(nat,indzmat,GeomTmp,GeomTmpC)

        GeomTmpC2=reshape(GeomTmpC,(/Nat,3/))

        DO I=1,Nat

           IF (renum) THEN

              Iat=Order(I)

              !             WRITE(IOCART,'(1X,A10,3(1X,F15.8))') Trim(AtName(I)),GeomTmpC(3*Iat-2:3*Iat)

              WRITE(IOCART,'(1X,A10,3(1X,F15.8))') Trim(AtName(I)),GeomTmpC2(Iat,1:3)

           ELSE

              Iat=OrderInv(I)

              !            WRITE(IOCART,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomTmpC(3*I-2:3*I)

              WRITE(IOCART,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomTmpC2(I,1:3)

           END IF

        END DO

     END IF ! matches IF (COORD.EQ.'ZMAT') THEN

     IF (COORD.EQ.'BAKER') THEN

        ! Writing the Cartesian coordinates

        WRITE(IOCART,'(1X,I5)') Nat

        WRITE(IOCART,'(1X,A)') Title

        ! we have to generate the cartesian coordinates from the internal coordinates

        ! we can also print the cartesiann coordinates directly from XyzGeomF??

        !Call ConvertBakerInternal_cart(InCoordF(IGeom-1,:),InCoordF(IGeom,:))

        ! XyzGeomF(NGeomF,3,Nat), GeomTmpC(3*Nat),GeomTmpC2(Nat,3)

        !GeomTmpC2=reshape(GeomTmpC,(/Nat,3/)). XyzGeomF is already updated in 

        ! EgradPath.f90

        DO I =1, 3

           GeomTmpC2(:,I)=XyzGeomF(IGeom,I,:)

        END DO

        DO I=1,Nat

           IF (renum) THEN

              Iat=Order(I)

              !             WRITE(IOCART,'(1X,A10,3(1X,F15.8))') Trim(AtName(I)),GeomTmpC(3*Iat-2:3*Iat)

              !Print *, GeomTmpC2(Iat,1:3)

              WRITE(IOCART,'(1X,A10,3(1X,F15.8))') Trim(AtName(I)),GeomTmpC2(Iat,1:3)

           ELSE

              Iat=OrderInv(I)

              !            WRITE(IOCART,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomTmpC(3*I-2:3*I)

              WRITE(IOCART,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomTmpC2(I,1:3)

           END IF

        END DO

     END IF ! matches IF (COORD.EQ.'BAKER') THEN

     IF (COORD.EQ.'MIXED') THEN

        ! Writing the Cartesian coordinates

        WRITE(IOCART,'(1X,I5)') Nat

        WRITE(IOCART,'(1X,A)') Title

        ! we have to generate the cartesian coordinates from the internal coordinates

        Call Mixed2Cart(nat,indzmat,geomtmp,GeomTmpC)

        GeomTmpC2=reshape(GeomTmpC,(/Nat,3/))

        DO I=1,Nat

           IF (renum) THEN

              Iat=Order(I)

              !             WRITE(IOCART,'(1X,A10,3(1X,F15.8))') Trim(AtName(I)),GeomTmpC(3*Iat-2:3*Iat)

              WRITE(IOCART,'(1X,A10,3(1X,F15.8))') Trim(AtName(I)),GeomTmpC2(Iat,1:3)

           ELSE

              Iat=OrderInv(I)

              !             WRITE(IOCART,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomTmpC(3*I-2:3*I)

              WRITE(IOCART,'(1X,A10,3(1X,F15.8),A)') Trim(AtName(Iat)),GeomTmpC2(I,1:3)

           END IF

        END DO

     END IF ! matches IF (COORD.EQ.'MIXED') THEN

  END DO ! matches DO IGeom=1,NGeomF

  CLOSE(IOTMP)

  IF ((COORD.EQ.'ZMAT').OR.(COORD.EQ.'MIXED').OR.(COORD.EQ.'BAKER')) CLOSE(IOCART)

  DEALLOCATE(GeomTmpC,GeomTmp)

  if (debug) Call Header('Exiting Write_PATH')

END SUBROUTINE Write_path