SUBROUTINE Calc_Zmat_frag(na,atome,ind_zmat,val_zmat,x,y,z,r_cov,fact)

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!   Goal: to compute a viable Z-Matrix starting from the 
!         cartesian coordinates of the atoms
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Input:
!  na    : Number or atoms 
! atome  : Mass number of the atoms (H=1, He=2,...)
! x,y,z  : cartesian coordinates of the atoms
! r_cov  : array containing the covalent radii of the atoms
!  fact  : multiplicative factor used to determine if two atoms are linked.
!          see CalcCnct for more details.
!
! Output:
! ind_zmat : INTEGER(na,5) contains the indices of the Z-matrix
! val_zmat : REAL(na,3)    contains the values of the Z-matrix
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! History
!
! v1.0 written for Cart a long time ago. Does not use fragment analysis, but was quite good !
!
! v2.0 12/2007
! Comes from Calc_zmat_constr_frag, based on a fragment analysis of the 
! system under study.
! It should be more flexible and robust.
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!----------------------------------------------------------------------
!  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, only : max_Z, NMaxL, Nom, Pi
  Use Io_module

  IMPLICIT NONE

  integer(KINT), INTENT(IN) ::  na,atome(na)
  real(KREAL), INTENT(IN)   ::  x(Na),y(Na),z(Na),fact
  real(KREAL), INTENT(IN)   ::  r_cov(0:Max_Z)

  INTEGER(KINT), INTENT(OUT) :: ind_zmat(Na,5)
  real(KREAL), INTENT(OUT)   :: val_zmat(Na,3)

  INTEGER(KINT) :: idx_zmat(NA)
! Frozen contains the indices of frozen atoms
! INTEGER(KINT) Frozen(*),NFroz
! Number of fragment, Index of the current fragment for loops
  INTEGER(KINT) :: NbFrag
! Fragment(I) contains the fragment index to which I belongs
! NbAtFrag(j) contains the number of atoms of fragment j
  INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !na
! FragAt(i,:) lists the atoms of fragment i
  INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !na,na
! MaxLFrag(i,1) contains the maximum of links for an atom for fragment i
! MaxLFrag(i,2) is the atom that has this number of linked atoms
  INTEGER(KINT), ALLOCATABLE :: MaxLFrag(:,:) !na,2
! INTEGER(KINT), ALLOCATABLE :: IdxFragAt(:) !na
! INTEGER(KINT), ALLOCATABLE :: FrozBlock(:,:) !(na,0:na)
! DistFrag contains the distance between a given atom and some other atoms
  REAL(KREAL), ALLOCATABLE :: DistFrag(:) !na
! FragDist(I) contains the index of the atoms corresponding to DistFrag(I)
  INTEGER(KINT), ALLOCATABLE :: FragDist(:) ! na

  INTEGER(KINT) :: IdxCaFaire, IAfaire
  INTEGER(KINT), ALLOCATABLE ::  LIAISONS(:,:) ! (Na,0:NMaxL)
! INTEGER(KINT), ALLOCATABLE ::  LiaisonsIni(:,:) ! (Na,0:NMaxL)
  INTEGER(KINT), ALLOCATABLE :: CAFaire(:) ! (Na)

  real(KREAL) ::  vx1,vy1,vz1,norm1
  real(KREAL) ::  vx2,vy2,vz2,norm2
  real(KREAL) ::  vx3,vy3,vz3,norm3
  real(KREAL) ::  vx4,vy4,vz4,norm4
  real(KREAL) ::  vx5,vy5,vz5,norm5
  real(KREAL) ::  val,val_d, d12, d13,d23,d
  Logical Debug, FirstAt, DebugGaussian
  LOGICAL, ALLOCATABLE :: DejaFait(:), FCaf(:) !(na)
!  Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen
  LOGICAL F1213, F1223,F1323

  INTEGER(KINT) :: I,J,n0,n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax
  INTEGER(KINT) :: I3, I1, Ip
  INTEGER(KINT) :: I0, Izm, JAt
  INTEGER(KINT) :: OrderZmat

  REAL(KREAL) :: sDihe
  REAL(KREAL), PARAMETER :: LocalNCart=0.d0

  INTERFACE
     function valid(string) result (isValid)
       CHARACTER(*), intent(in) :: string
       logical                  :: isValid
     END function VALID

     FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
       real(KREAL) ::  v1x,v1y,v1z,norm1
       real(KREAL) ::  v2x,v2y,v2z,norm2
       real(KREAL) ::  angle
     END FUNCTION ANGLE

     FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
       real(KREAL) ::  v1x,v1y,v1z,norm1
       real(KREAL) ::  v2x,v2y,v2z,norm2
       real(KREAL) ::  SinAngle
     END FUNCTION SINANGLE

     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
       real(KREAL) ::  v1x,v1y,v1z,norm1
       real(KREAL) ::  v2x,v2y,v2z,norm2
       real(KREAL) ::  v3x,v3y,v3z,norm3
       real(KREAL) ::  angle_d,ca,sa
     END FUNCTION ANGLE_D
  END INTERFACE

  debug=valid("calc_zmat").OR.valid("calc_zmat_frag")
  debugGaussian=valid("zmat_gaussian")
  
  if (debug) WRITE(*,*) "=============================== Entering Calc_zmat_frag ========================"
  if (na.le.2) THEN
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
     RETURN
  END IF

  ALLOCATE(FragDist(Na),Fragment(na), NbAtFrag(na),FragAt(na,na))
! ALLOCATE(FrozFrag(na,3), IdxFragAt(na), FrozBlock(na,0:na))
! ALLOCATE(FrozFrag(na,3))
  ALLOCATE(DistFrag(na),Liaisons(na,0:NMaxL))
! ALLOCATE(Liaisons(na,0:NMaxL),LiaisonsIni(na,0:NMaxL))
! ALLOCATE(CaFaire(na),DejaFait(Na),FCaf(na),FrozAt(na))
  ALLOCATE(CaFaire(na+1),DejaFait(Na),FCaf(na))

  if (debug) THEN
     WRITE(*,*) "DBG Calc_zmat_frag - Cartesian coordinates"
     DO I=1,na
        WRITE(*,'(1X,A3,3(1X,F15.8))') Nom(atome(i)),x(i),y(i),z(i)
     END DO
  END if

  DO I=1,na
     DO J=1,5
        Ind_Zmat(I,J)=0
     END DO
  END DO
  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!     Easy case : 3 or less atoms
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  if (Na.eq.3) THEN
     d12=sqrt((x(1)-x(2))**2+(y(1)-y(2))**2+(z(1)-z(2))**2)
     d13=sqrt((x(1)-x(3))**2+(y(1)-y(3))**2+(z(1)-z(3))**2)
     d23=sqrt((x(3)-x(2))**2+(y(3)-y(2))**2+(z(3)-z(2))**2)
     F1213=(d12<=d13)
     F1323=(d13<=d23)
     F1223=(d12<=d23)
     if (debug) THEN
        WRITE(*,*) "DEBUG Calc_Zmat 3 atoms"
        WRITE(*,*) "d12,d13,d23:",d12,d13,d23
        WRITE(*,*) "F1213,F1323,F1223:",F1213,F1323,F1223
     END IF
     OrderZmat=0
     if (F1213) orderZmat=OrderZmat+4
     if (F1323) orderZmat=OrderZmat+2
     if (F1223) orderZmat=OrderZmat+1
     if (debug) WRITE(*,*) 'OrderZmat=',OrderZmat
     SELECT CASE (OrderZmat)
        CASE (0)
! F F F ordre 2-3----1 
           ind_zmat(1,1)=3
           ind_zmat(2,1)=2
           ind_zmat(2,2)=3
           ind_zmat(3,1)=1
           ind_zmat(3,2)=3
           ind_zmat(3,3)=2
        CASE (2)
! F T F ordre 1-3----2
           ind_zmat(1,1)=3
           ind_zmat(2,1)=1
           ind_zmat(2,2)=3
           ind_zmat(3,1)=2
           ind_zmat(3,2)=3
           ind_zmat(3,3)=1
        CASE (3)
! F T T ordre 2---1-3
           ind_zmat(1,1)=1
           ind_zmat(2,1)=3
           ind_zmat(2,2)=1
           ind_zmat(3,1)=2
           ind_zmat(3,2)=1
           ind_zmat(3,3)=3
        CASE (5)
! T F T  ordre 1-2----3
           ind_zmat(1,1)=2
           ind_zmat(2,1)=1
           ind_zmat(2,2)=2
           ind_zmat(3,1)=3
           ind_zmat(3,2)=2
           ind_zmat(3,3)=1
        CASE (7)
! T T T ordre 3----1-2
           ind_zmat(1,1)=1
           ind_zmat(2,1)=2
           ind_zmat(2,2)=1
           ind_zmat(3,1)=3
           ind_zmat(3,2)=1
           ind_zmat(3,3)=2
        END SELECT

     IF (debug) THEN
        WRITE(*,*) "DBG Calc_zmat_frag - Nat=3 -"
        DO i=1,na
           WRITE(*,'(1X,A3,5(1X,I5))') Nom(Atome(ind_zmat(i,1))),(ind_zmat(i,j),j=1,5)
        END DO
     END IF

     ! We have ind_zmat, we fill val_zmat
     val_zmat=0.d0
     n2=ind_zmat(2,1)
     n1=ind_zmat(2,2)
     d=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2+(z(n1)-z(n2))**2)
     val_zmat(2,1)=d
     n1=ind_zmat(3,1)
     n2=ind_zmat(3,2)
     n3=ind_zmat(3,3)
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
     if (debug) write(*,*) n1,n2,norm1
     val_zmat(3,1)=norm1

     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
     if (debug) write(*,*) n2,n3,norm2,val
     val_zmat(3,2)=val

     RETURN
  END IF !matches if (Na.eq.3) THEN
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  End of Easy case : 3 or less atoms
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! General Case 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!

! Initialization
  DejaFait=.False.
  Liaisons=0
  ind_zmat=0
  val_zmat=0.d0

  if (debug) WRITE(*,*) "Coucou from Calc_zmat_frag.f90; L240"
  
  if (debug) THEN
     WRITE(*,*) "Bonds initialized"
     WRITE(*,*) 'Covalent radii used'
     DO iat=1,na
        i=atome(iat)
        WRITE(*,*) Nom(I),Iat,r_cov(i),r_cov(i)*fact
     END DO
  END IF

1003 FORMAT(1X,I4,' - ',25(I5))

!   First step : connectivity  are calculated

  Call CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)

  if (debug) THEN
     WRITE(*,*) "Connections calculated"
     DO IL=1,na
        WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
     END DO
  END IF

  FCaf=.TRUE.
  Call Decomp_frag(na,liaisons,FCaf,nbfrag,Fragment,NbAtFrag,FragAt)
  
  IF (debug) THEN
     WRITE(*,*) 'I found ',NbFrag, 'fragments'
     WRITE(*,*) (NbAtFrag(I),I=1,NbFrag)
     DO I=1,NbFrag
        WRITE(*,*) NbAtFrag(I)
        WRITE(*,*) 'Fragment ', i
        DO J=1,Na
           IF (Fragment(J).EQ.I) WRITE(*,'(1X,I3,1X,A5,3(1X,F10.6))') J,Nom(Atome(J)) &
                ,X(J),Y(J),Z(J)
        END DO
        WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
     END DO
  END IF

 ALLOCATE(MaxLFrag(NbFrag,2))

 MaxLFrag=0

  DO I=1,NbFrag
     MaxLFrag(I,1)=Liaisons(FragAt(I,1),0)
     MaxLFrag(I,2)=FragAt(I,1)

     DO J=1,NbAtFrag(I)
      Iat=FragAt(I,J)
        IF (Liaisons(IAt,0).GT.MaxLFrag(I,1)) THEN
       MaxLFrag(I,1)=Liaisons(Iat,0)
       MaxLFrag(I,2)=Iat
    END IF
     END DO
     IF (debug) WRITE(*,*) 'Frag :',I,', atom ',MaxLFrag(I,2), ' has ',MaxLFrag(I,2),' links'
  END DO

  !     We will now build the molecule fragment by fragment
  !     We choose the starting fragment with two criteria:
  !     1- Number of linked atoms:
  !       * >=3 is good as it fully defines the coordinate space
  !       * 2 is ok as we can either use a 3rd atom from the same fragment
  !       or add a X atom somewhere but this complicates quite a lot the way
  !       to treat the conversion from cartesian to zmat latter
  !       * 1 is bad...
  !     2- Size of the fragment
  !     this allows us to deal more easily with cases 1- when number of 
  !     directly linked atoms is less than 3

  IFrag=0
! I0 is the number of connections of the best fragment
  I0=0
! I1 is the number of atoms of the best fragment
  I1=0
  IAt=0
  DO I=1,NbFrag
    SELECT CASE(MaxLFrag(I,1)-I0)
     CASE (1:)
        IFrag=I
        I0=MaxLFrag(I,1)
        I1=NbAtFrag(I)
        IAt=MaxLFrag(I,2)
     CASE (0)
        if (NbAtFrag(I).GT.I1) THEN
           IFrag=I
           I0=MaxLFrag(I,1)
           I1=NbAtFrag(I)
           IAt=MaxLFrag(I,2)
        END IF
     END SELECT
  END DO

  if (debug) WRITE(*,'(1X,A,I5,A,I5,A,I5,A,I5)') 'Starting with fragment:',IFrag,' with ',I0   &
      ,' max links for atom',IAt,' fragment size',NbAtFrag(IFrag)

  !     We will build the first fragment in a special way, as it will
  !     set the coordinates system

  if (debug) WRITE(*,*) 'Fragment 1, starting with atom:',IAt, &   
    'with ',I0,' connections'

  DejaFait=.FALSE.
  FCaf=.FALSE.
  
  izm=0
  SELECT CASE (I0)
  CASE(3:)
     if (debug) WRITE(*,*) 'DBG select case I0 3'
     n0=Iat

     ITmp=2
     sDihe=0.
     n2=IAt
     n3=Liaisons(Iat,1)
     ! We search for the third atom while making sure that it is not aligned with first two !
     DO While ((ITmp.LE.Liaisons(Iat,0)).AND.(sDihe.LE.0.09d0))
        n4=Liaisons(Iat,Itmp)
        CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
        CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
        val_d=angle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
        sDiHe=abs(sin(val_d*pi/180.d0))
        if (debug) Write(*,*) 'Trying n2,n3,n4,sdihe,val_d',n2,n3,n4,sdihe,val_d
        Itmp=Itmp+1
     END DO
     If (debug) WRITE(*,*) 'Itmp,Liaisons',Itmp,Liaisons(Iat,1:NMaxL)
     Liaisons(Iat,Itmp-1)=Liaisons(iat,2)
     Liaisons(Iat,2)=n4
     If (debug) WRITE(*,*) 'Itmp,Liaisons',Itmp,Liaisons(Iat,1:NMaxL)

     if (sDihe.LE.0.09d0) THEN
        WRITE(*,*) "PROBLEM !!! All atoms linked to ",n0," are aligned..."
    WRITE(*,*) "Surprising as this atom has at least three bonds... For NOW STOP"
    STOP
     END IF

     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,  vx5,vy5,vz5,norm5)


! We search for the fourth atom while checking that it is not aligned with 1st and 2nd atoms.
     Itmp=3
     sDiHe=0.
! PFL 28 Dec 2007 ->
! I had a test on the dihedral angle, but I cannot see why it is important to have
! a non planar fragment at the begining... ethylene works and is fully planar
! I thus suppress this test
!
!     DO While ((ITmp.LE.Liaisons(Iat,0)).AND.(sDihe.LE.0.09d0))
!        ITmp=ITmp+1
        n1=Liaisons(Iat,Itmp)
        if (debug) WRITe(*,*) 'trying n1,n2,n3,n4',n1,n2,n3,n4
        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
        ! Is this atom aligned with n2-n3 ?
        val_d=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
        sDiHe=abs(sin(val_d*pi/180.d0))
        if (debug) Write(*,*) 'Angle n3-n2-n1',val_d
        if (sDiHe.le.0.09d0) THEN
           ! As atoms n2,n3 and n4 are not aligned, we interchange n3 and n4 so that n1,n2 and n3 are not aligned
           if (debug) WRITE(*,*) "n3-n2-n1 aligned, we interchange n3 and n4"
           n1=n3
           n3=n4
           n4=n1
           n1=Liaisons(Iat,ITmp)
           CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)           
           CALL vecteur(n2,n4,x,y,z,vx3,vy3,vz3,norm3)           
           val_d=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
           if (debug) Write(*,*) 'NEW Angle n3-n2-n1',val_d
        END IF

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        ! avant de tester l'angle diedre, il faut verifier que ce 4e atome n'est pas
        ! aligne avec les 2 premiers. 
        ! comme on a deja testé que les 3 premiers ne sont pas alignes,
        ! s'il est aligné avec les 2 premiers, on peut inverser le role de 2 et 3.
        ! On pourrait tout simplier en faisant une bete recherche parmi tous les atomes geles
        ! de ce bloc (au moins 4 ?) avec le critere 1) on les range par distance croissante
        ! 2) on les scanne tant que l'angle valence n'est pas bon, puis tant que diehedre pas bon
        ! on pourrait comme ca faire un tableau avec les atomes ranges d'abord pour le fragment sélectionné
        ! puis les atomes des autres fragment par distance croissante.
        ! les autres fragments ne seraient additonnés que si l'on ne trouve pas notre bonheur dans le premier bloc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
        val_d=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
             vx2,vy2,vz2,norm2)
        sDihe=abs(sin(val_d*pi/180.d0))
        if (debug) WRITE(*,*) 'n2,n3,n4,n1, angle_d',n2,n3,n4,n1,val_d
!     END DO

     DejaFait(n2)=.TRUE.
     DejaFait(n3)=.TRUE.
     DejaFait(n4)=.TRUE.

!     if (sDihe.LE.0.09d0) THEN
!        !     None of the atoms linked to IAt are good to define the third
!        !     direction in space...
!        !     We will look at the other atoms
!        !     we might improve the search so as to take the atom closest to IAt
!        if (debug) WRITE(*,*) "All atoms linked to ",Iat," are in a plane. Looking for other atoms"
!        ITmp=0
!        DO I=1,NbAtFrag(IFrag)
!           JAt=FragAt(IFrag,I)
!           if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
!              n1=JAt
!              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
!              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
!              val_d=angle_d(vx4,vy4,vz4,norm4, &
!                   vx5,vy5,vz5,norm5, &
!                   vx2,vy2,vz2,norm2)
!              if (abs(sin(val_d)).GE.0.09D0) THEN
!                 ITmp=ITmp+1
!                 DistFrag(ITmp)=Norm1
!                 FragDist(ITmp)=JAt
!              END IF
!           END IF
!        END DO

!        IF (ITmp.EQ.0) THEN
!           !     All dihedral angles between frozen atoms are less than 5°
!           !     (or more than 175°). We have to look at other fragments :-(
!           DO I=1,NFroz
!              Jat=Frozen(I)
!              if (.NOT.DejaFait(JAt)) THEN
!                 n1=JAt
!                 CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
!                 CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
!                 val_d=angle_d(vx4,vy4,vz4,norm4, &
!                      vx5,vy5,vz5,norm5, &
!                      vx2,vy2,vz2,norm2)
!                 if (abs(sin(val_d)).GE.0.09D0) THEN
!                    ITmp=ITmp+1
!                    DistFrag(ITmp)=Norm1
!                    FragDist(ITmp)=JAt
!                 END IF
!              END IF
!           END DO
!           IF (ITmp.EQ.0) THEN
!              !     All frozen atoms are in a plane... too bad
!              WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
!              WRITE(*,*) 'For now, I do not treat this case'
!              STOP
!           END IF
!        END IF
!        I1=0
!        d=1e9
!        DO I=1,ITmp
!           IF (DistFrag(I).LE.d) THEN
!              d=DistFrag(I)
!              I1=FragDist(I)
!           END IF
!        END DO
!     ELSE                !(sDihe.LE.0.09d0)
!        I1=FrozBlock(IAt,ITmp)
!        if (debug) WRITE(*,*) 'I1,n1:',I1,n1
!     END IF                 !(sDihe.LE.0.09d0)
!     !     we now have the atom that is closer to the first one and that
!     !     forms a non 0/Pi dihedral angle
!
!  <- PFL 28 Dec 2007

! We construct the begining of the Z-Matrix

     ind_zmat(1,1)=n2
     ind_zmat(2,1)=n3
     ind_zmat(2,2)=n2
     ind_zmat(3,1)=n4
     ind_zmat(3,2)=n2
     ind_zmat(3,3)=n3
     DejaFait(n2)=.TRUE.
     DejaFait(n3)=.TRUE.
     DejaFait(n4)=.TRUE.
     CaFaire(1)=n3
     CaFaire(2)=n4

! PFL 28 Dec 2007
! We have not selected a fourth atom, so that the following is not needed
!     ind_zmat(4,1)=I1
!     ind_zmat(4,2)=n2
!     ind_zmat(4,3)=n3
!     ind_zmat(4,4)=n4
!     DejaFait(I1)=.TRUE.
!     CaFaire(3)=I1
!     CaFaire(4)=0
!     IdxCaFaire=4
!     izm=4
!     FCaf(I1)=.TRUE.
!!!!!!!
! and replaced by:
     CaFaire(3)=0
   IdxCaFaire=3
   izm=3
!
! <- PFL 28 Dec 2007
   
     FCaf(n2)=.TRUE.
     FCaf(n3)=.TRUE.
   FirstAt=.TRUE.
     DO I=3,Liaisons(Iat,0)
        IF (.NOT.DejaFait(Liaisons(Iat,I))) THEN
           izm=izm+1
           !           ind_zmat(izm,1)=Liaisons(Iat,I)
           !           ind_zmat(izm,2)=n2
           !           ind_zmat(izm,3)=n3
           !           ind_zmat(izm,4)=n4
           Call add2indzmat(na,izm,Liaisons(Iat,I),n2,n3,n4,ind_zmat,x,y,z)
           if (FirstAt) THEN
           n4=Liaisons(Iat,I)
         FirstAt=.FALSE.
       END IF
           IF (.NOT.FCaf(Liaisons(Iat,I))) THEN
              CaFaire(IdxCaFaire)=Liaisons(Iat,I)
              IdxCaFaire=IdxCaFaire+1
              CaFaire(IdxCaFaire)=0
              FCaf(Liaisons(Iat,I))=.TRUE.
           END IF
           DejaFait(Liaisons(Iat,I))=.TRUE.
        END IF
     END DO

     if (debug) THEN
        WRITE(*,*) "Ind_zmat 0 - SELECT CASE I0 3: -- izm=",izm
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
        DO I=4,izm
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2),  &
                ind_zmat(I,3), ind_zmat(I,4)
        END DO
     END IF


     !     First four atoms (at least) have been put we go on with next parts
     !     of this fragment... later


  CASE(2)
     if (debug) WRITE(*,*) 'DBG select case I0 2'
   WRITE(*,*) "PFL 28 Dec 2007: No test of alignment here :(  TO DO TO DO TO DO"
     ind_zmat(1,1)=IAt
     ind_zmat(2,1)=Liaisons(IAt,1)
     ind_zmat(2,2)=IAt
     ind_zmat(3,1)=Liaisons(IAt,2)
     ind_zmat(3,2)=IAt
     ind_zmat(3,3)=Liaisons(IAt,1)
     DejaFait(IAt)=.TRUE.
     DejaFait(Liaisons(Iat,1))=.TRUE.
     DejaFait(Liaisons(Iat,2))=.TRUE.
     CaFaire(1)=Liaisons(Iat,1)
     CaFaire(2)=Liaisons(Iat,2)
     FCaf(Liaisons(Iat,1))=.TRUE.
     FCaf(Liaisons(Iat,2))=.TRUE.

! PFL 28 Dec 2007 ->
! We do NOT need a fourth atom !!! The third direction in space is defined by the cross
! product of the first two directions
!
! the following is thus commented
!
!     !     We search for a fourth atom, first in the FrozBlock atoms
!     ITmp=2
!     sDihe=0.
!     n2=IAt
!     n3=Liaisons(Iat,1)
!     n4=Liaisons(Iat,2)
!     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
!     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
!     CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2, vx5,vy5,vz5,norm5)
!
!     DO While ((ITmp.LE.FrozBlock(Iat,0)).AND.(sDihe.LE.0.09d0))
!        ITmp=ITmp+1
!        n1=FrozBlock(Iat,Itmp)
!        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
!        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
!        val_d=angle_d(vx4,vy4,vz4,norm4,  &
!             vx5,vy5,vz5,norm5,           &
!             vx2,vy2,vz2,norm2)
!        sDihe=abs(sin(val_d))
!     END DO
!     IF (debug) WRITE(*,*) 'DBG 4th atom, ITmp, sDihe',ITmp, sDihe
!     if (sDihe.LE.0.09d0) THEN
!        !     None of the frozen atoms linked to IAt are good to define the third
!        !     direction in space...
!        !     We will look at the other frozen atoms
!        !     we might improve the search so as to take the atom closest to IAt
!        ITmp=0
!        DO I=1,NbAtFrag(IFrag)
!           JAt=FragAt(IFrag,I)
!           if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
!              n1=JAt
!              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
!              CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
!              val_d=angle_d(vx4,vy4,vz4,norm4,    &
!                   vx5,vy5,vz5,norm5,              &
!                   vx2,vy2,vz2,norm2)
!              if (abs(sin(val_d)).GE.0.09D0) THEN
!                 ITmp=ITmp+1
!                 DistFrag(ITmp)=Norm1
!                 FragDist(ITmp)=JAt
!              END IF
!           END IF
!        END DO
!        IF (ITmp.EQ.0) THEN
!           !     All dihedral angles between frozen atoms are less than 5°
!           !     (or more than 175°). We have to look at other fragments :-(
!           DO I=1,NFroz
!              Jat=Frozen(I)
!              if (.NOT.DejaFait(JAt)) THEN
!                 n1=JAt
!                 CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
!                 CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
!                 val_d=angle_d(vx4,vy4,vz4,norm4,       &
!                      vx5,vy5,vz5,norm5,                 &
!                      vx2,vy2,vz2,norm2)
!                 if (abs(sin(val_d)).GE.0.09D0) THEN
!                    ITmp=ITmp+1
!                    DistFrag(ITmp)=Norm1
!                    FragDist(ITmp)=JAt
!                 END IF
!              END IF
!           END DO
!           IF (ITmp.EQ.0) THEN
!              !     All frozen atoms are in a plane... too bad
!              WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
!              WRITE(*,*) 'For now, I do not treat this case'
!              STOP
!           END IF
!        END IF
!        I1=0
!        d=1e9
!        DO I=1,ITmp
!           IF (DistFrag(I).LE.d) THEN
!              d=DistFrag(I)
!              I1=FragDist(I)
!           END IF
!        END DO
!     ELSE                   !(sDihe.LE.0.09d0)
!        I1=FrozBlock(IAt,ITmp)
!     END IF                 !(sDihe.LE.0.09d0)
!     !     we now have the atom that is closer to the first one and that
!     !     forms a non 0/Pi dihedral angle
!     !     ind_zmat(4,1)=I1
!     !     ind_zmat(4,2)=IAt
!     !     ind_zmat(4,3)=Liaisons(Iat,1)
!     !     ind_zmat(4,4)=Liaisons(Iat,2)
!     n3=Liaisons(Iat,1)
!     n4=Liaisons(Iat,2)
!     Call add2indzmat(na,4,I1,Iat,n3,n4,ind_zmat,x,y,z)
!     Liaisons(Iat,1)=n3
!     Liaisons(Iat,2)=n4
!     DejaFait(I1)=.TRUE.
!     CaFaire(3)=I1
!     CaFaire(4)=0
!     IdxCaFaire=4
!     izm=4
!     FCaf(I1)=.TRUE.
!
!!!!!! <- PFL 28 Dec 2007

     CaFaire(3)=0
     IdxCaFaire=3
     izm=3

  CASE(1)
     if (debug) WRITE(*,*) 'DBG select case I0 1, NbAtFrag=',NbAtFrag(IFrag)
     ind_zmat(1,1)=IAt
     ind_zmat(2,1)=Liaisons(IAt,1)
     ind_zmat(2,2)=IAt
     DejaFait(IAt)=.TRUE.
     DejaFait(Liaisons(Iat,1))=.TRUE.
     IdxCaFaire=2
     CaFaire(1)=Liaisons(Iat,1)
     CaFaire(2)=0
     FCaf(Liaisons(Iat,1))=.TRUE.

! PFL 28 Dec 2007 ->
! We do NOT need a fourth atom. So we will look only for a third atom
!
!!!!
!
     !     We search for a third and fourth atoms, first in the FrozBlock atoms
     !     It should not be possible to have   (FrozBlock(Iat,0).GT.2) and
     !     iat linked to only one atom !       


     !     we calculate the distances between Iat and all other frozen
     !     atoms of this fragment, and store only those for which 
     !     valence angles are not too close to 0/Pi. (limit:5°)

     ITmp=0
     CALL vecteur(Liaisons(Iat,1),IAt,x,y,z,vx2,vy2,vz2,norm2)

! PFL 28 Dec 2007: As MaxL=1 I think that there is at most 2 atoms in this fragment... 
! so that the following loop is useless... this should be tested more carefully
     DO I=1,NbAtFrag(IFrag)
        JAt=FragAt(IFrag,I)
        if (.NOT.DejaFait(JAt)) THEN
           CALL vecteur(JAt,IAt,x,y,z,vx1,vy1,vz1,norm1)
           if (abs(cos(angle(vx1,vy1,vz1,norm1,  &
                vx2,vy2,vz2,norm2))).LE.0.996d0) THEN
              ITmp=ITmp+1
              DistFrag(ITmp)=Norm1
              FragDist(ITmp)=JAt
           END IF
        END IF
     END DO

     IF (ITMP.EQ.0) THEN
        !     We have no atoms correct in this fragment, we use
        !     atoms from other fragments
        DO Jat=1,Na
! DejaFait(Iat)=.TRUE. so that we do not need to test Jat/=Iat
           if (.NOT.DejaFait(JAt)) THEN
              CALL vecteur(JAt,IAt,x,y,z,vx1,vy1,vz1,norm1)
              if (abs(cos(angle(vx1,vy1,vz1,norm1, &
                   vx2,vy2,vz2,norm2))).LE.0.996d0) THEN
                 ITmp=ITmp+1
                 DistFrag(ITmp)=Norm1
                 FragDist(ITmp)=JAt
              END IF
           END IF
        END DO
        IF (ITMP.EQ.0) THEN
           WRITE(*,*) 'It seems all atoms are aligned'
           WRITE(*,*) 'Case non treated for now :-( '
           STOP
        END IF
     END IF

     I1=0
     d=1e9
! PFL 28 Dec 2007: There exists some F90 intrinsics to find the smallest element of an array.
! The following loop should be replaced by it !
     DO I=1,ITmp
        IF (DistFrag(I).LE.d) THEN
           I1=FragDist(I)
           d=DistFrag(I)
        END IF
     END DO

     !     we now have the atom that is closer to the first one and that
     !     forms a non 0/Pi valence angle
     ind_zmat(3,1)=I1
     ind_zmat(3,2)=IAt
     ind_zmat(3,3)=Liaisons(Iat,1)
     DejaFait(I1)=.TRUE.
     CaFaire(2)=I1
     FCaf(I1)=.TRUE.


! PFL 28 Dec 2007 ->
! We do NOT need a fourth atom so that the following is suppressed
!
!     !     we search for a fourth atom !
!     !     We search for a fourth atom, first in the FrozBlock atoms
!     ITmp=2
!     sDihe=0.
!     n2=IAt
!     n3=Liaisons(Iat,1)
!     n4=I1
!     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
!     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
!     CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2,vx5,vy5,vz5,norm5)
!
!     !     We will look at the other frozen atoms
!     !     we might improve the search so as to take the atom closest to IAt
!     ITmp=0
!     DO I=1,NbAtFrag(IFrag)
!        JAt=FragAt(IFrag,I)
!        if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
!           n1=JAt
!           CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
!           CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
!           val_d=angle_d(vx4,vy4,vz4,norm4,  &
!                vx5,vy5,vz5,norm5,   &
!                vx2,vy2,vz2,norm2)
!           if (abs(sin(val_d)).GE.0.09D0) THEN
!              ITmp=ITmp+1
!              DistFrag(ITmp)=Norm1
!              FragDist(ITmp)=JAt
!           END IF
!        END IF
!     END DO
!     IF (ITmp.EQ.0) THEN
!        !     All dihedral angles between frozen atoms are less than 5°
!        !     (or more than 175°). We have to look at other fragments :-(
!        DO I=1,NFroz
!           Jat=Frozen(I)
!           if (.NOT.DejaFait(JAt)) THEN
!              n1=JAt
!              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
!              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
!              val_d=angle_d(vx4,vy4,vz4,norm4,   &
!                   vx5,vy5,vz5,norm5,           &
!                   vx2,vy2,vz2,norm2)
!              if (abs(sin(val_d)).GE.0.09D0) THEN
!                 ITmp=ITmp+1
!                 DistFrag(ITmp)=Norm1
!                 FragDist(ITmp)=JAt
!              END IF
!           END IF
!       END DO
!        IF (ITmp.EQ.0) THEN
!           !     All frozen atoms are in a plane... too bad
!           WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
!           WRITE(*,*) 'For now, I do not treat this case'
!           STOP
!        END IF
!     END IF                 ! ITmp.EQ.0 after scaning fragment
!     I1=0
!     d=1e9
!     DO I=1,ITmp
!        IF (DistFrag(I).LE.d) THEN
!           d=DistFrag(I)
!           I1=FragDist(I)
!        END IF
!     END DO
!
!     !     we now have the atom that is closer to the first one and that
!     !     forms a non 0/Pi dihedral angle
!     !     ind_zmat(4,1)=I1
!     !     ind_zmat(4,2)=IAt
!     !     ind_zmat(4,3)=ind_zmat(2,1)
!     !     ind_zmat(4,4)=ind_zmat(3,1)
!     n3=ind_zmat(2,1)
!     n4=ind_zmat(3,1)
!     Call add2indzmat(na,4,I1,IAt,n3,n4,ind_zmat,x,y,z)
!     ind_zmat(2,1)=n3
!     ind_zmat(3,1)=n4
!     DejaFait(I1)=.TRUE.
!     CaFaire(3)=I1
!     CaFaire(4)=0
!     IdxCaFaire=4
!     izm=4
!     FCaf(I1)=.TRUE.
!!!!!!!!!!!
!
! <- PFL 28 Dec 2007

     CaFaire(3)=0
   IdxCaFaire=3
   
  CASE(0)
     WRITE(*,*) 'All atoms are separated .. '
     WRITE(*,*) 'this case should be treated separately !'
     STOP
  END SELECT

  if (debug) THEN
     WRITE(*,*) 'ind_zmat 1 izm=',izm
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
     DO I=4,izm
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2),  &
             ind_zmat(I,3), ind_zmat(I,4)
     END DO
  END IF

  DO I=1,izm
     Idx_zmat(ind_zmat(I,1))=i
  END Do

  !     at least first three atoms of this fragment done...
  !     we empty the 'cafaire' array before going on
  IAFaire=1
  DO WHILE (CaFaire(IaFaire).NE.0)
     n1=CaFaire(IaFaire)
     n2=ind_zmat(idx_zmat(N1),2)
     if (idx_zmat(N1).EQ.2) THEN
        !     We have a (small) problem: we have to add atoms linked to the 2nd
        !     atom of the zmat. This is a pb because we do not know
        !     which atom to use to define the dihedral angle
        !     we take the third atom of the zmat
        n3=ind_zmat(3,1)
     ELSE
        n3=ind_zmat(idx_zmat(n1),3)
     END IF
   
   FirstAt=.TRUE.
     DO I=1,Liaisons(n1,0)
        IAt=Liaisons(n1,I)
! PFL 29.Aug.2008 ->
! We dissociate the test on 'DejaFait' that indicates that this atom
! has already been put in the Zmat
! from the test on FCaf that indicates that this atom has been put in the
! 'CAFaire' list that deals with identifying its connectivity.
! Those two test might differ in some cases.
        IF (.NOT.DejaFait(Iat)) THEN
           izm=izm+1
           if (debug) WRITE(*,*) ">1< Adding atom ",Iat,"position izm=",izm
           !           ind_zmat(izm,1)=iat
           !           ind_zmat(izm,2)=n1
           !           ind_zmat(izm,3)=n2
           !           ind_zmat(izm,4)=n3
           Call add2indzmat(na,izm,iat,n1,n2,n3,ind_zmat,x,y,z)
       if (FirstAt) THEN
          n3=Iat
        FirstAt=.FALSE.
           END IF
           idx_zmat(iat)=izm
           DejaFait(iat)=.TRUE.
        END IF
        IF (.NOT.FCaf(Iat)) THEN
           CaFaire(IdxCaFaire)=iat
           IdxCaFaire=IdxCaFaire+1
           CaFaire(IdxCaFaire)=0
           FCaf(Iat)=.TRUE.
        END IF
! <- PFL 29.Aug.2008
     END DO
     IaFaire=IaFaire+1
  END Do                    ! DO WHILE CaFaire

  if (debug) THEN
     WRITE(*,*) 'ind_zmat 2, izm=',izm
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
     DO I=4,izm
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2),  &
             ind_zmat(I,3), ind_zmat(I,4)
     END DO
  END IF

  !     We have finished putting atoms linked to the first one
  ! There should not be any atom left from this fragment. We check:
  !     we will add other atoms of this fragment
  DO I=1,NbAtFrag(IFrag)
     Iat=FragAt(IFrag,I)
     if (debug) WRITE(*,*) "DBG: I,Iat,dejafait",I,Iat,DejaFait(Iat)
     IF (.NOT.DejaFait(Iat)) THEN
           WRITE(*,*) 'Treating atom I,Iat',I,Iat
        
     END IF

  END DO

  NbAtFrag(Ifrag)=0
  MaxLFrag(IFrag,1)=0
  
  !     we start again with the rest of the molecule...
  ! v 1.01 We add the fragment in the order corresponding to NbAtFrag
  KMax=NbFrag-1

  IF (DEBUG) WRITE(*,*) "Adding the ",Kmax," remaining fragments"
  DO K=1, KMax
  IFrag=0
  I0=0
  IAt=0
  I1=0
  DO I=1,NbFrag
    SELECT CASE(MaxLFrag(I,1)-I0)
     CASE (1:)
        IFrag=I
        I0=MaxLFrag(I,1)
        IAt=MaxLFrag(I,2)
        I1=NbAtFrag(I)
     CASE (0)
        if (NbAtFrag(I).GT.I1) THEN
           IFrag=I
           I0=MaxLFrag(I,1)
           IAt=MaxLFrag(I,2)
           I1=NbAtFrag(I)
        END IF
     END SELECT
 
  END DO

  if (debug) WRITE(*,'(1X,A,I5,A,I5,A,I5,A,I5)') 'Adding fragment:',IFrag,' with ',I0   &
      ,' max links for atom',IAt,' fragment size',NbAtFrag(IFrag)
    
     MaxLFrag(IFrag,1)=0

! We search for the closest atoms of the previous fragments to the atom with max links
  d=1e9
  DO J=1,izm
    Call vecteur(iat,ind_zmat(j,1),x,y,z,vx1,vy1,vz1,norm1)
    if (norm1.le.d) THEN
      Jat=j
      d=norm1
    END IF
  END DO
  n2=ind_zmat(jat,1)
  SELECT CASE(jat)
        CASE (1)
     n3=ind_zmat(2,1)
     n4=ind_zmat(3,1)
    CASE (2)
     n3=ind_zmat(1,1)
     n4=ind_zmat(3,1)
    CASE DEFAULT
     n3=ind_zmat(jAt,2)
     n4=ind_zmat(jat,3)
  END SELECT
  izm=izm+1
  Call add2indzmat(na,izm,iat,n2,n3,n4,ind_zmat,x,y,z)
    idx_zmat(iat)=izm
    DejaFait(iat)=.TRUE.
  IdxCaFaire=2
    CaFaire(1)=iat
              CaFaire(2)=0
              FCaf(Iat)=.TRUE.
              IaFaire=1
              DO WHILE (CaFaire(IaFaire).NE.0)
                 n1=CaFaire(IaFaire)
                 n2=ind_zmat(idx_zmat(N1),2)
                 if (idx_zmat(N1).EQ.2) THEN
                    !     We have a (small) problem: we have to add atoms linked to the 2nd
                    !     atom of the zmat. This is a pb because we do not know
                    !     which atom to use to define the dihedral angle
                    !     we take the third atom of the zmat
                    n3=ind_zmat(3,1)
                 ELSE
                    n3=ind_zmat(idx_zmat(n1),3)
                 END IF
                 DO I3=1,Liaisons(n1,0)
                    IAt=Liaisons(n1,I3)
! PFL 29.Aug.2008 ->
! We dissociate the test on 'DejaFait' that indicates that this atom
! has already been put in the Zmat
! from the test on FCaf that indicates that this atom has been put in the
! 'CAFaire' list that deals with identifying its connectivity.
! Those two test might differ for a frozen atom linked to non frozen atoms.
                    IF (.NOT.DejaFait(Iat)) THEN
                       izm=izm+1
                       Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
                       idx_zmat(iat)=izm
                       n3=iat
                       DejaFait(Iat)=.TRUE.
                    END IF
                    IF (.NOT.FCaf(Iat)) THEN
                       CaFaire(IdxCaFaire)=iat
                       IdxCaFaire=IdxCaFaire+1
                       CaFaire(IdxCaFaire)=0
                       FCaf(Iat)=.TRUE.
                    END IF
! <- PFL 29.Aug.2008
                 END DO
                 IaFaire=IaFaire+1
              END Do              ! DO WHILE CaFaire

        if (debug) THEN
           WRITE(*,*) 'ind_zmat 4'
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
           DO Ip=4,izm
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1), ind_zmat(Ip,2), &
                   ind_zmat(Ip,3), ind_zmat(Ip,4)
           END DO
        END IF

  END DO                    ! Loop on all fragments


     ! We have ind_zmat. We calculate val_zmat :-)
     if (debug) WRITE(*,*) "Calculating val_zmat"

     val_zmat(1,1)=0.d0
     val_zmat(1,2)=0.d0
     val_zmat(1,3)=0.d0
     val_zmat(2,2)=0.d0
     val_zmat(2,3)=0.d0
     val_zmat(3,3)=0.d0

     n1=ind_zmat(2,1)
     n2=ind_zmat(2,2)

     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)

     val_zmat(2,1)=norm1


     n1=ind_zmat(3,1)
     n2=ind_zmat(3,2)
     n3=ind_zmat(3,3)

     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)

     val_zmat(3,1)=norm1

     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)

     val_zmat(3,2)=val

     DO i=4,na

        n1=ind_zmat(i,1)
        n2=ind_zmat(i,2)
        n3=ind_zmat(i,3)
        n4=ind_zmat(i,4)

        if (debug) WRITE(*,*) "Doing i,n1,n2,n3,n4",i,n1,n2,n3,n4
        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)

        CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
        val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)

        CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
        CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)

        val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &
             vx2,vy2,vz2,norm2)

        !     write(*,11) n1,n2,norm1,n3,val,n4,val_d
        !11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))

        val_zmat(i,1)=norm1
        val_zmat(i,2)=val
        val_zmat(i,3)=val_d

     END DO

     if (debug) THEN
        WRITE(*,*) 'DBG Cre_Zmat_Frag: ind_zmat'
        DO I=1,na
           WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
        END DO

        WRITE(*,*) 'DBG Cre_Zmat_Frag: Full zmat'
        DO I=1,na
           WRITE(*,'(1X,I5,1X,I5,F8.4,2(1X,I5,1X,F7.2))') ind_zmat(i,1),(ind_zmat(i,j+1),val_zmat(i,j),j=1,3)
        END DO

     END IF

     if (debugGaussian) THEN
        WRITE(*,*) 'DBG Cre_Zmat_Frag: Gaussian Zmat - START'
        Call WriteMixed_Gaussian(na,atome,LocalNCart,ind_zmat,val_zmat)
        WRITE(*,*) 'DBG Cre_Zmat_Frag: Gaussian Zmat - END'
     END IF


     if (debug) WRITE(*,*) "Deallocate (FragDist,Fragment, NbAtFrag,FragAt)"
     DEALLOCATE(FragDist,Fragment, NbAtFrag,FragAt)
     if (debug) WRITE(*,*) "Deallocate (DistFrag,Liaisons)"
     DEALLOCATE(DistFrag,Liaisons)
     if (debug) WRITE(*,*) "Deallocate(CaFaire,DejaFait)"
     DEALLOCATE(CaFaire,DejaFait,FCaf,MaxLFrag)
 


     if (debug) WRITE(*,*) "=============================== Exiting Calc_zmat_frag ========================"

END SUBROUTINE Calc_Zmat_frag

