SUBROUTINE Calc_Zmat_const_frag(na,atome,ind_zmat,val_zmat,x,y,z, &
     r_cov,fact,frozen)

  !     
  !     This second version enables the user to freeze some atoms
  !     the frozen atoms indices are listed in the frozen array.
  !     
  !     The idea is surely too stupid to work all the time... but here it is
  !     we first construct the zmat using only the frozen atoms, and then
  !     we add the other atoms on top of the first ones...
  !     
  Use Path_module, only : max_Z, NMaxL, Nom,MaxFroz
  Use Io_module

  IMPLICIT NONE


  CHARACTER(5) :: AtName
  integer(KINT) :: na,atome(na),ind_zmat(Na,5)
  INTEGER(KINT) :: idx_zmat(NA)
  real(KREAL) ::  x(Na),y(Na),z(Na),fact
  real(KREAL) ::  val_zmat(Na,3)
  real(KREAL) :: r_cov(0:Max_Z)
  !     Frozen contains the indices of frozen atoms
  INTEGER(KINT) :: Frozen(*),NFroz
  INTEGER(KINT), ALLOCATABLE :: FrozDist(:) ! na
  INTEGER(KINT) :: NbFrag,IdxFrag
  INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !na
  INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !na,na
  INTEGER(KINT), ALLOCATABLE :: FrozFrag(:,:) !na,3
  !  INTEGER(KINT), ALLOCATABLE :: IdxFragAt(:) !na
  INTEGER(KINT), ALLOCATABLE :: FrozBlock(:,:) !(na,0:na)
  REAL(KREAL), ALLOCATABLE :: DistFroz(:) !na

  INTEGER(KINT) :: IdxCaFaire, IAfaire
  INTEGER(KINT), ALLOCATABLE ::  LIAISONS(:,:) ! (Na,0:NMaxL)
  INTEGER(KINT), ALLOCATABLE ::  LiaisonsBis(:,:) ! (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 PasFini,Debug
  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) :: IMax, I3,I1,Ip, IFragFroz
  INTEGER(KINT) :: I0, Izm, JAt,Izm0

  REAL(KREAL) :: sDihe, Pi

  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)
       use Path_module, only :  Pi,KINT, KREAL
       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)
       use Path_module, only :  Pi,KINT, KREAL
       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)
       use Path_module, only :  Pi,KINT, KREAL
       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



  Pi=dacos(-1.d0)
  debug=valid("calc_zmat_constr").OR.valid("calc_zmat_contr_frag")

  if (debug) WRITE(*,*) "=============================== Entering Calc_zmat_contr_frag ========================"
  if (na.le.2) THEN
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
     RETURN
  END IF

  ALLOCATE(FrozDist(Na),Fragment(na), NbAtFrag(na),FragAt(na,na))
  !  ALLOCATE(FrozFrag(na,3), IdxFragAt(na), FrozBlock(na,0:na))
  ALLOCATE(FrozFrag(na,3), FrozBlock(na,0:na))
  ALLOCATE(DistFroz(na),Liaisons(na,0:NMaxL))
  ALLOCATE(LiaisonsBis(na,0:NMaxL),LiaisonsIni(na,0:NMaxL))
  ALLOCATE(CaFaire(na),DejaFait(Na),FCaf(na),FrozAt(na))

  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
     if (F1213) THEN
        if (F1323) THEN
           ind_zmat(1,1)=3
           ind_zmat(2,1)=1
           ind_zmat(2,2)=3
           ind_zmat(3,1)=2
           ind_zmat(3,2)=1
           ind_zmat(3,3)=3
        ELSE
           ind_zmat(1,1)=1
           ind_zmat(2,1)=2
           ind_zmat(2,2)=1
           ind_zmat(3,1)=3
           ind_zmat(3,2)=2
           ind_zmat(3,3)=1
        END IF
     ELSE
        IF (F1223) THEN
           ind_zmat(1,1)=2
           ind_zmat(2,1)=1
           ind_zmat(2,2)=2
           ind_zmat(3,1)=3
           ind_zmat(3,2)=1
           ind_zmat(3,3)=2
        ELSE
           ind_zmat(1,1)=2
           ind_zmat(2,1)=3
           ind_zmat(2,2)=2
           ind_zmat(3,1)=1
           ind_zmat(3,2)=3
           ind_zmat(3,3)=2
        END IF
     END IF
     IF (debug) THEN
        DO i=1,na
           WRITE(*,'(1X,5(1X,I5))') (ind_zmat(i,j),j=1,5)
        END DO
     END IF

     !     We have ind_zmat, we fill val_zmat
     val_zmat(1,1)=0.
     val_zmat(1,2)=0.
     val_zmat(1,3)=0.
     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
     val_zmat(2,2)=0.
     val_zmat(2,3)=0.
     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
     val_zmat(3,3)=0.

     RETURN
  END IF
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!  End of   Easy case : 3 or less atoms
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! General Case 
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! 1 - Frozen Atoms

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

  NFroz=0
  I=1
  DO WHILE (Frozen(I).NE.0)
     If (Frozen(I).LT.0) THEN
        DO J=Frozen(I-1),abs(Frozen(I))
           IF (.NOT.FrozAt(J)) THEN
              NFroz=NFroz+1
              FrozAt(J)=.TRUE.
           END IF
        END DO
     ELSEIF (.NOT.FrozAt(Frozen(I))) THEN
        FrozAt(Frozen(I))=.TRUE.
        NFroz=NFroz+1
     END IF
  I=I+1
  END DO

  if (debug) WRITE(*,*) 'DGB zmtc NFroz=', NFroz
  if (debug) WRITE(*,*) 'DGB zmtc FrozAt=',(FrozAt(I),I=1,na)

  if (debug) THEN
     WRITE(*,*) "Liaisons 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))
  !     DO IL=1,na
  !     WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
  !     WRITE(IOOUT,1003) Il,(LIAISONSBis(IL,JL),JL=0,NMaxL)
  !     END DO
!   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(*,'(20(1X,I4))') (NbAtFrag(I),I=1,NbFrag)
     DO I=1,NbFrag
        WRITE(*,*) NbAtFrag(I)
        WRITE(*,*) 'Fragment ', i
        DO J=1,Na
           IF (Fragment(J).EQ.I)   THEN                
            if (FrozAt(J))    THEN   
                   WRITE(*,'(1X,I3,"f",1X,A5,3(1X,F10.6))') J,Nom(Atome(J)) &
                        ,X(J),Y(J),Z(J)
               ELSE
                    WRITE(*,'(1X,I3,2X,A5,3(1X,F10.6))') J,Nom(Atome(J)) &
                      ,X(J),Y(J),Z(J)
                END IF
       END IF
        END DO
        WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
     END DO
  END IF


! First, we decompose the connectivity into Frozen atoms and non frozen atoms
  DO I=1,na
     LiaisonsBis(I,0)=0
     DO J=1,Liaisons(I,0)
        IF (FrozAt(Liaisons(I,J))) THEN
           LiaisonsBis(I,0)=LiaisonsBis(I,0)+1
           LiaisonsBis(I,LiaisonsBis(I,0))=Liaisons(I,J)
        END IF
     END DO
     IMax=LiaisonsBis(I,0)+1
     LiaisonsBis(I,Imax)=0
     DO J=1,Liaisons(I,0)
        IF (.NOT.FrozAt(Liaisons(I,J))) THEN
           LiaisonsBis(I,IMax)=LiaisonsBis(I,Imax)+1
           LiaisonsBis(I,LiaisonsBis(I,Imax)+Imax)=Liaisons(I,J)
        END IF
     END DO
  END DO

  if (debug) THEN
     WRITE(*,*) "Liaisons and LiaisonsBis"
     DO I=1,Na
        WRITE(*,'(1X,I3," -",I3,":",15(1X,I3))') I, &
             (Liaisons(I,J),J=0,NMaxL)
        WRITE(*,'(1X,I3," +",I3,":",15(1X,I3))') I, &
             (LiaisonsBis(I,J),J=0,NMaxL)
     END DO
  END IF

! Now, for each frozen atom, we count the length of the frozen block
! FrozBlock(I,0) contains the number of frozen atoms forming a frozen
!                fragment containing atom I
! FrozBlock(I,:) is the list of the frozen atoms of this fragment.
  Do I=1,na
     FrozBlock(I,0)=-1
  END DO
  DO I=1,na
     IF (FrozAt(I).AND.(FrozBlock(I,0).EQ.-1)) THEN
        if (debug) WRITE(*,*) 'Treating atom ',I
        FrozBlock(I,0)=1
        FrozBlock(I,1)=I
        DO J=1,na
           DejaFait(J)=.FALSE.
        END DO
        DejaFait(I)=.TRUE.
        DO J=1,LiaisonsBis(I,0)
           CaFaire(J)=LiaisonsBis(I,J)
        END DO
        IdxCaFaire=LiaisonsBis(I,0)+1
        CaFaire(IdxCaFaire)=0
        IAfaire=1
        FCaf=DejaFait
        DO WHILE (CaFaire(IAfaire).NE.0)
           Iat=CaFaire(IAFAire)
           if (debug) WRITE(*,*) 'IaFaire; Iat:', IaFaire, Iat, DejaFait(Iat)
           IF (.NOT.DejaFait(Iat)) THEN
              FrozBlock(I,0)=FrozBlock(I,0)+1
              FrozBlock(I,FrozBlock(I,0))=Iat
              DO J=1,LiaisonsBis(Iat,0)
                 IF ((.NOT.DejaFait(LiaisonsBis(Iat,J))).AND.(.NOT.FCaf(LiaisonsBis(Iat,J)))) THEN
                    CaFaire(IdxCaFaire)=LiaisonsBis(Iat,J)
                    IdxCaFaire=IdxCaFaire+1
                    CaFaire(IdxCaFaire)=0
                    FCaf(LiaisonsBis(Iat,J))=.TRUE.
                 END IF
              END DO
              !     WRITE(*,*) 'liaisonbis(Iat,0),FrozBlick(I,0)',&
              !                      LiaisonsBis(Iat,0), FrozBlock(I,0)
           END IF
           DejaFait(Iat)=.TRUE.
           IaFaire=IaFaire+1
        END DO
        if (debug) WRITE(*,*) 'I,FrozBlock(0),IaFaire',I,FrozBlock(I,0),IaFaire
        if (debug) WRITE(*,*) 'FrozBlock:',FrozBlock(I,1:FrozBlock(I,0))
        !        FrozBlock(I,1)=I
        !        DO J=2,IaFaire
        !           FrozBlock(I,J)=CaFaire(J-1)
        !        END DO
        !     We copy this block into all frozen atoms that belongs to it
        DO J=2,Frozblock(I,0)
           Iat=FrozBlock(I,J)
           DO K=0,FrozBlock(I,0)
              FrozBlock(Iat,K)=FrozBlock(I,K)
           END DO
        END DO
     ELSE
        IF (.NOT.FrozAt(I))      FrozBlock(I,0)=0
     END IF
  END DO

  if (debug) THEN
     WRITE(*,*) "FrozBlock"
     DO I=1,Na
        If (FrozAt(I))  WRITE(*,'(1X,I3," -",I3,":",15(1X,I3))') I, &
             (FrozBlock(I,J),J=0,FrozBlock(I,0))
     END DO
  END IF

  DO I=1,NbFrag
     FrozFrag(I,1)=0
     FrozFrag(I,2)=0
     DO J=1,NbAtFrag(I)
        IF (FrozAt(FragAt(I,J))) THEN
           FrozFrag(I,1)=FrozFrag(I,1)+1
           IF (FrozBlock(FragAt(I,J),0).GE.FrozFrag(I,2))  &
                FrozFrag(I,2)=FrozBlock(FragAt(I,J),0)
           if (FrozFrag(I,3).LE.LiaisonsBis(FragAt(I,J),0))& 
                FrozFrag(I,3)=LiaisonsBis(FragAt(I,J),0)
        END IF
     END DO
     IF (debug) WRITE(*,*) 'Frag :',I,FrozFrag(I,1), &
          ' frozen atoms,max linked:'  &
          ,FrozFrag(I,2),' with at most',FrozFrag(I,3),' frozen'
  END DO

  !     We will now build the molecule fragment by fragment
  !     First the frozen atoms, then the rest of the fragment
  !     We choose the starting fragment with two criteria:
  !     1- Number of linked frozen 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- Max number of linked frozen atoms
  !     this allows us to deal more easily with cases 1- when number of 
  !     directly linked frozen atoms is less than 3

  IFrag=0
  I0=0
  I1=0
  DO I=1,NbFrag
     SELECT CASE(FrozFrag(I,3)-I0)
     CASE (1:)
        IFrag=I
        I0=FrozFrag(I,3)
        I1=FrozFrag(I,2)
     CASE (0)
        if (FrozFrag(I,2).GT.I1) THEN
           IFrag=I
           I0=FrozFrag(I,3)
           I1=FrozFrag(I,2)
        END IF
     END SELECT
  END DO

  if (debug) WRITE(*,'(1X,A,I5,A,I5,A,I5,A)') 'Starting with fragment:',IFrag,' with ',I0,' frozen linked and overall',I1,' linked'

  !     We will build the first fragment in a special way, as it will
  !     set the coordinates system
  !     We look for the frozen atom that is linked to the maximum frozen atom,
  !     and belongs to the longer fragment
  I0=0
  I1=0
  DO I=1,NbAtFrag(IFrag)
     IF (FrozAt(FragAt(IFrag,I))) THEN
        SELECT CASE(LiaisonsBis(FragAt(IFrag,I),0)-I0)
        CASE (1:) 
           IAt=FragAt(IFrag,I)
           I0=LiaisonsBis(IAt,0)
           I1=FrozBlock(IAt,0)
        CASE (0)
           IF (FrozBlock(FragAt(IFrag,I),0).GT.I1) THEN
              IAt=FragAt(IFrag,I)
              I0=LiaisonsBis(Iat,0)
              I1=FrozBlock(Iat,0)
           END IF
        END SELECT
     END IF
     if (debug) WRITE(*,*) 'DBG: I,FragAt(IFrag,I),IAt,I0,I1',I,FragAt(IFrag,I),IAt,I0,I1
  END DO

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

  DO I=1,na
     DejaFait(I)=.FALSE.
     FCaf(I)=.FALSE.
  END DO

  izm=0
  SELECT CASE (I0)
  CASE(3:)
     if (debug) WRITE(*,*) 'DBG select case I0 3'
     n0=Iat
     !     We search for the fourth atom while making sure that the dihedral angle
     !     is not 0 or Pi
     ITmp=2
     sDihe=0.
     n2=IAt
     n3=LiaisonsBis(Iat,1)
     ! We search for the third atom while making sure that it is not aligned with first two !
     DO While ((ITmp.LE.LiaisonsBis(Iat,0)).AND.(sDihe.LE.0.09d0))
        n4=LiaisonsBis(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,LiaisonsBis',Itmp,LiaisonsBis(Iat,1:NMaxL)
     LiaisonsBis(Iat,Itmp-1)=LiaisonsBis(iat,2)
     LiaisonsBis(Iat,2)=n4
     If (debug) WRITE(*,*) 'Itmp,LiaisonsBis',Itmp,LiaisonsBis(Iat,1:NMaxL)

     if (sDihe.LE.0.09d0) THEN
        WRITE(*,*) "Dans le caca car tous les atoms de ce block sont alignés!"
     END IF

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


     Itmp=2
     sDiHe=0.

     DO While ((ITmp.LT.LiaisonsBis(Iat,0)).AND.(sDihe.LE.0.09d0))
        ITmp=ITmp+1
        n1=LiaisonsBis(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=LiaisonsBis(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,norm1,vx2,vy2,vz2,norm2, &
             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 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
        if (debug) WRITE(*,*) "All atoms linked to ",Iat," are in a plane. Looking for other frozen 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,norm1,vx2,vy2,vz2,norm2, &
                   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
                 DistFroz(ITmp)=Norm1
                 FrozDist(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,norm1, &
                      vx2,vy2,vz2,norm2, &
                      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
                    DistFroz(ITmp)=Norm1
                    FrozDist(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 (DistFroz(I).LE.d) THEN
              d=DistFroz(I)
              I1=FrozDist(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

     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

     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

     FCaf(n2)=.TRUE.
     FCaf(n3)=.TRUE.
     FCaf(I1)=.TRUE.
     izm=4
     DO I=3,LiaisonsBis(Iat,0)
        IF (.NOT.DejaFait(LiaisonsBis(Iat,I))) THEN
           izm=izm+1
           !           ind_zmat(izm,1)=LiaisonsBis(Iat,I)
           !           ind_zmat(izm,2)=n2
           !           ind_zmat(izm,3)=n3
           !           ind_zmat(izm,4)=n4
           Call add2indzmat(na,izm,LiaisonsBis(Iat,I),n2,n3,n4,ind_zmat,x,y,z)
           IF (.NOT.FCaf(LiaisonsBis(Iat,I))) THEN
              CaFaire(IdxCaFaire)=LiaisonsBis(Iat,I)
              IdxCaFaire=IdxCaFaire+1
              CaFaire(IdxCaFaire)=0
              FCaf(LiaisonsBis(Iat,I))=.TRUE.
           END IF
           DejaFait(LiaisonsBis(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'
     if (debug) WRITE(*,*) 'ATTENTION For now, case with only 3 frozen atoms non treated'
     ind_zmat(1,1)=IAt
     ind_zmat(2,1)=liaisonsBis(IAt,1)
     ind_zmat(2,2)=IAt
     ind_zmat(3,1)=LiaisonsBis(IAt,2)
     ind_zmat(3,2)=IAt
     ind_zmat(3,3)=LiaisonsBis(IAt,1)
     DejaFait(IAt)=.TRUE.
     DejaFait(LiaisonsBis(Iat,1))=.TRUE.
     DejaFait(LiaisonsBis(Iat,2))=.TRUE.
     CaFaire(1)=LiaisonsBis(Iat,1)
     CaFaire(2)=LiaisonsBis(Iat,2)
     FCaf(LiaisonsBis(Iat,1))=.TRUE.
     FCaf(LiaisonsBis(Iat,2))=.TRUE.

     !     We search for a fourth atom, first in the FrozBlock atoms
     ITmp=2
     sDihe=0.
     n2=IAt
     n3=LiaisonsBis(Iat,1)
     n4=LiaisonsBis(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,norm3, &
          vx2,vy2,vz2,norm2, &
          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,norm1, &
             vx2,vy2,vz2,norm2, &
             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,norm1,  &
                   vx2,vy2,vz2,norm2,               &
                   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
                 DistFroz(ITmp)=Norm1
                 FrozDist(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,norm1,   &
                      vx2,vy2,vz2,norm2,                 &
                      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
                    DistFroz(ITmp)=Norm1
                    FrozDist(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 (DistFroz(I).LE.d) THEN
              d=DistFroz(I)
              I1=FrozDist(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)=LiaisonsBis(Iat,1)
     !     ind_zmat(4,4)=LiaisonsBis(Iat,2)
     n3=LiaisonsBis(Iat,1)
     n4=LiaisonsBis(Iat,2)
     Call add2indzmat(na,4,I1,Iat,n3,n4,ind_zmat,x,y,z)
     LiaisonsBis(Iat,1)=n3
     LiaisonsBis(Iat,2)=n4
     DejaFait(I1)=.TRUE.
     CaFaire(3)=I1
     CaFaire(4)=0
     IdxCaFaire=4
     izm=4
     FCaf(I1)=.TRUE.

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

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

     IF (FrozBlock(Iat,0).GT.2) THEN
        WRITE(*,*) 'I found some inconsistancy: IAt linked to 1'
        WRITE(*,*) 'Atom only, but belongs to a frozblock of at '
        WRITE(*,*) 'least 3 atoms. IAt, FrozBlock'
        WRITE(*,*) Iat,(FrozBlock(IAt,J),J=0,FrozBlock(Iat,0))
        STOP
     END IF

     !     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(LiaisonsBis(Iat,1),IAt,x,y,z,vx2,vy2,vz2,norm2)

     DO I=1,NbAtFrag(IFrag)
        JAt=FragAt(IFrag,I)
        if (FrozAt(Jat).AND.(.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
              DistFroz(ITmp)=Norm1
              FrozDist(ITmp)=JAt
           END IF
        END IF
     END DO

     IF (ITMP.EQ.0) THEN
        !     We have no frozen atoms correct in this fragment, we use
        !     atoms from other fragments
        DO I=1,NFroz
           Jat=Frozen(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
                 DistFroz(ITmp)=Norm1
                 FrozDist(ITmp)=JAt
              END IF
           END IF
        END DO
        IF (ITMP.EQ.0) THEN
           WRITE(*,*) 'It seems all frozen atoms are aligned'
           WRITE(*,*) 'Case non treated for now :-( '
           STOP
        END IF
     END IF

     I1=0
     d=1e9
     DO I=1,ITmp
        IF (DistFroz(I).LE.d) THEN
           I1=FrozDist(I)
           d=DistFroz(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)=LiaisonsBis(Iat,1)
     DejaFait(I1)=.TRUE.
     CaFaire(2)=I1
     FCaf(I1)=.TRUE.


     !     we search for a fourth atom !
     !     We search for a fourth atom, first in the FrozBlock atoms
     ITmp=2
     sDihe=0.
     n2=IAt
     n3=LiaisonsBis(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,norm3,  &
          vx2,vy2,vz2,norm2,    &
          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,norm1, &
                vx2,vy2,vz2,norm2, &
                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
              DistFroz(ITmp)=Norm1
              FrozDist(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,norm1,  &
                   vx2,vy2,vz2,norm2, &
                   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
                 DistFroz(ITmp)=Norm1
                 FrozDist(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 (DistFroz(I).LE.d) THEN
           d=DistFroz(I)
           I1=FrozDist(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.

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

  if (debug) THEN
     WRITE(*,*) 'ind_zmat 1'
     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 four (frozen) atoms of this fragment done...
  !     we empty the 'cafaire' array before pursuing
  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
     IF (LiaisonsBis(n1,0).GE.1) THEN
        IAt=LiaisonsBis(n1,1)
        if (.NOT.DejaFait(Iat)) THEN
           izm=izm+1
           if (debug) WRITE(*,*) ">0< 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)
           Idx_zmat(iat)=izm
           n3=iat
           IF (.NOT.FCaf(Iat)) THEN
              CaFaire(IdxCaFaire)=iat
              IdxCaFaire=IdxCaFaire+1
              CaFaire(IdxCaFaire)=0
              FCaf(Iat)=.TRUE.
           END IF
           DejaFait(iat)=.TRUE.
        END IF
     END IF
     DO I=2,LiaisonsBis(n1,0)
        IAt=LiaisonsBis(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 for a frozen atom linked to non frozen atoms.
        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)
           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
  !     we will add other frozen atoms of this fragment
  DO I=1,NbAtFrag(IFrag)
     Iat=FragAt(IFrag,I)
     if (debug) WRITE(*,*) "DBG: I,Iat,Frozat,dejafait,frozbloc",I,Iat,FrozAt(Iat), DejaFait(Iat),FrozBlock(Iat,0)
     IF (FrozAt(Iat).AND.(.NOT.DejaFait(Iat))) THEN
        !     in order to have the zmat as connected as possible,
        !     we look if this atom belongs to a frozblock
        if (debug) WRITE(*,*) 'Treating atom I,Iat, FrozBlock ',I,Iat, FrozBlock(Iat,0)
        IF (FrozBlock(Iat,0).EQ.1) THEN
           !     it is a lonely atom :-)
           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
           !           ind_zmat(izm,1)=iat
           !           ind_zmat(izm,2)=n2
           !           ind_zmat(izm,3)=n3
           !           ind_zmat(izm,4)=n4
           Call add2indzmat(na,izm,iat,n2,n3,n4,ind_zmat,x,y,z)
           DejaFait(iat)=.TRUE.
           idx_zmat(iat)=iat
        ELSE        !(FrozBlock(Iat,0).EQ.1)
           !     we have a group of atoms
           !     We search for the atom of the group with the most links
           ITmp=-1
           DO J=1,FrozBlock(Iat,0)
              Jat=FrozBlock(Iat,J)
              IF ((.NOT.DejaFait(Jat)).AND.  &
                   (LiaisonsBis(Jat,0).GT.ITMP)) THEN
                 JL=Jat
                 ITmp=LiaisonsBis(Jat,0)
              END IF
           END DO
           if (debug) WRITE(*,*) 'JL,ITmp:',JL,ITmp
           Iat=JL
           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
           if (debug) WRITE(*,*) 'Jat,d:',Jat,d
           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
           !           ind_zmat(izm,1)=iat
           !           ind_zmat(izm,2)=n2
           !           ind_zmat(izm,3)=n3
           !           ind_zmat(izm,4)=n4
           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.

           if (debug) WRITE(*,*) izm,iat,n2,n3,n4

           IaFaire=1
           DO WHILE (CaFaire(IaFaire).NE.0)
              n1=CaFaire(IaFaire)
              n2=ind_zmat(idx_zmat(N1),2)
              if (debug) WRITE(*,*) 'DBG Cafaire, Iafaire,n1,n2',Iafaire,n1,n2
              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
              if (debug) WRITe(*,*) 'DBG :n3,liaisonBis',n3,LiaisonsBis(n1,0)
              DO I3=1,LiaisonsBis(n1,0)
! 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.
                 IAt=LiaisonsBis(n1,I3)
                 if (.NOT.DejaFait(Iat)) THEN
                    izm=izm+1
                    !                 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)
                    idx_zmat(iat)=izm
                    if (I3.EQ.1) n3=ind_zmat(izm,1)
                    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
        END IF
     END IF
     if (debug) THEN
        WRITE(*,*) 'ind_zmat 3'
        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

  FrozFrag(IFrag,1)=-1
  FrozFrag(IFrag,2)=-1
  FrozFrag(IFrag,3)=-1

  !     we start again with the rest of the molecule...
  ! v 1.01 We add the fragment in the order corresponding to FrozFrag
  KMax=0
  DO I=1,NbFrag
     IF (FrozFrag(I,1).NE.0) KMax=KMax+1
  END DO

  IF (DEBUG) WRITE(*,*) "Adding the",Kmax,"fragments with frozen atoms"
  DO K=1, KMax-1
     IFrag=0
     I0=0
     I1=0
     DO I=1,NbFrag
        SELECT CASE(FrozFrag(I,3)-I0)
        CASE (1:)
           IFrag=I
           I0=FrozFrag(I,3)
           I1=FrozFrag(I,2)
        CASE (0)
           if (FrozFrag(I,2).GT.I1) THEN
              IFrag=I
              I0=FrozFrag(I,3)
              I1=FrozFrag(I,2)
           END IF
        END SELECT
     END DO

     FrozFrag(IFrag,1)=-1
     FrozFrag(IFrag,2)=-1
     FrozFrag(IFrag,3)=-1

     if (debug) WRITE(*,*) "Adding Fragment ",IFrag,K

     DO I=1,NbAtFrag(IFrag)
        Iat=FragAt(IFrag,I)
        IF (FrozAt(Iat).AND.(.NOT.DejaFait(Iat))) THEN
           !     in order to have the zmat as connected as possible,
           !     we look if this atom belongs to a frozblock
           IF (FrozBlock(Iat,0).EQ.1) THEN
              !     it is a lonely atom :-)
              if (debug) write(*,*) "DBG 4: Lonely atom",Iat
              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
              !              ind_zmat(izm,1)=iat
              !              ind_zmat(izm,2)=n2
              !              ind_zmat(izm,3)=n3
              !              ind_zmat(izm,4)=n4
              Call add2indzmat(na,izm,iat,n2,n3,n4,ind_zmat,x,y,z)
              idx_zmat(iat)=izm
              DejaFait(iat)=.TRUE.
           ELSE
              !     we have a group of atoms
              !     We search for the atom of the group with the most links
              if (debug) write(*,*) "DBG 4b: ",Iat," belong to frozblock",FrozBlock(Iat,0)
              ITmp=-1
              DO J=1,FrozBlock(Iat,0)
                 Jat=FrozBlock(Iat,J)
                 IF ((.NOT.DejaFait(Jat)).AND.         &
                      (LiaisonsBis(Jat,0).GT.ITMP)) THEN
                    JL=Jat
                    ITmp=LiaisonsBis(Jat,0)
                 END IF
              END DO
              Iat=JL
              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
              !              ind_zmat(izm,1)=iat
              !              ind_zmat(izm,2)=n2
              !              ind_zmat(izm,3)=n3
              !              ind_zmat(izm,4)=n4
              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,LiaisonsBis(n1,0)
                    IAt=LiaisonsBis(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
           END IF
        END IF

        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 atoms of fragment IFrag
  END DO                    ! Loop on all fragments

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!        General case
!
! 2 - we have all frozen atoms done... now comes the non frozen atoms
! and we should fragment the fragments !
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

! First, we get rid of bonds between frozen atoms
! the trick is to do this on liaisons while keeping LiaisonsBis ok.

  if (debug) THEN
     WRITE(*,*) 'Frozen',NFroz
     WRITE(*,'(20(1X,I3))') (Frozen(I),I=1,NFroz)
  END IF

  DO I=1,na
     DO J=0,NMAxL
        LiaisonsIni(I,J)=LiaisonsBis(I,J)
     END DO
! PFL 29.Aug.2008 -> 
! We re-initialize FCaf in order to add frozen atoms that are connected
! to non frozen atoms
     FCaf(I)=.FALSE.
! <- 29.Aug.2008
  END DO

  DO I=1,Na
     IF (FrozAt(I)) THEN
        Iat=I
        if (debug) WRITE(*,'(20(1X,I3))') I,Iat,  &
             (LiaisonsIni(Iat,Itmp),ITmp=0,LiaisonsIni(Iat,0))
        DO J=1,LiaisonsIni(Iat,0)
           Jat=LiaisonsIni(Iat,J)
           Call Annul(Liaisons,Iat,Jat)
           Call Annul(Liaisons,Jat,Iat)
           Call Annul(LiaisonsIni,Jat,Iat)
           Liaisons(Iat,0)=Liaisons(Iat,0)-1
           Liaisons(Jat,0)=Liaisons(Jat,0)-1
           LiaisonsIni(Jat,0)=LiaisonsIni(Jat,0)-1
        END DO
     END IF
  END DO

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



  ! We re-decompose the system into fragments, but we omit on purpuse
  ! fragments consisting only of frozen atoms
  ! Now, frozblock will contain the frozen atom indices in a given fragment

  DO I=1,na
     Fragment(I)=0 
     NbAtFrag(I)=0
     FrozBlock(I,0)=0
  END DO
  IdxFrag=0
  NbFrag=0

  DO I=1,na
     IdxCAfaire=1
     CaFaire(IdxCaFaire)=0

     if (debug) WRITE(*,*) 'Treating atom I, fragment(I)',I,Fragment(I)
     IF (.NOT.FrozAt(I).OR.(Liaisons(I,0).NE.0)) THEN 
        IF (Fragment(I).EQ.0) THEN
           IdxFrag=IdxFrag+1
           NbFrag=NbFrag+1
           IFrag=IdxFrag
           Fragment(I)=IFrag
           NbAtFrag(IFrag)=NbAtFrag(IFrag)+1
           FragAt(IFrag,NbAtFrag(IFrag))=I
        ELSE  ! fragment(I).EQ.0
           IFrag=Fragment(I)
        END IF     ! fragment(I).EQ.0
        DO J=1,Liaisons(I,0)
           Iat=Liaisons(I,J)
           IF ((Fragment(Iat).NE.0).AND.(Fragment(Iat).NE.IFrag)) THEN
              WRITE(*,*) 'Error : Atoms ',I,' and ',Iat
              WRITE(*,*) 'are linked, but belongs to fragments ',IFrag,' and ',Fragment(Iat)
              STOP
           END IF
           IF (Fragment(Iat).EQ.0) THEN
              IF (.NOT.FCaf(IAt)) THEN
                 CaFaire(IdxCaFaire)=Iat
                 IdxCAfaire=IdxCAFaire+1
                 FCaf(IAt)=.TRUE.
              END IF
              Fragment(Iat)=IFrag
              NbAtFrag(IFrag)=NbAtFrag(IFrag)+1
              FragAt(IFrag,NbAtFrag(IFrag))=Iat
           END IF
        END DO
        CaFaire(IdxCaFaire)=0

        If (debug) WRITE(*,'(1X,A,1000I4)') 'Cafaire:',CaFaire(1:IdxCaFaire)
        If (debug) WRITE(*,*) 'IFrag=',IFrag

        IAfaire=1
        DO WHILE (CaFaire(IAfaire).NE.0)
           Iat=CaFaire(IaFaire)
           IF (.NOT.FrozAt(Iat).OR.(Liaisons(Iat,0).NE.0)) THEN 
              if (debug) WRITE(*,*) 'Cafaire treating ',Iat
              IF (Fragment(Iat).EQ.0) THEN
                 WRITE(*,*) 'Error: Atom ',Iat,' does not belong to any fragment !'
                 STOP
              END IF

              DO J=1,Liaisons(Iat,0)
                 Jat=Liaisons(Iat,J)
                 IF ((Fragment(Jat).NE.0).AND.(Fragment(Jat).NE.IFrag)) THEN
                    WRITE(*,*) 'Error: Atoms ',Iat,' and ',Liaisons(Iat,J)
                    WRITE(*,*) 'are linked, but belongs to fragments ',IFrag,' and ',Fragment(Liaisons(Iat,J))
                    STOP
                 END IF
                 IF (Fragment(Jat).EQ.0) THEN
                    IF (.NOT.FCaf(Liaisons(Iat,J))) THEN
                       CaFaire(IdxCaFaire)=Liaisons(Iat,J)
                       IdxCAfaire=IdxCAFaire+1
                       FCaf(Liaisons(Iat,J))=.TRUE.
                    END IF
                    Fragment(Jat)=IFrag
                    NbAtFrag(IFrag)=NbAtFrag(IFrag)+1
                    FragAt(IFrag,NbAtFrag(IFrag))=Jat
                 END IF
              END DO
              CaFaire(IdxCaFaire)=0
              If (debug) WRITE(*,*) 'IAfaire, IdxCaFaire,Cafaire :',IAfaire,IdxCafaire,' == ',CaFaire(IaFaire+1:IdxCaFaire)
              IAFaire=IAFaire+1
           END IF
        END DO
     END IF
  END DO

  ! We compute FrozBlock now
  DO IFrag=1,NbFrag
     DO I=1,NbAtFrag(IFrag)
        Iat=FragAt(IFrag,I)
        IF (FrozAt(Iat)) THEN
           FrozBlock(IFrag,0)=FrozBlock(IFrag,0)+1
           FrozBlock(IFrag,FrozBlock(IFrag,0))=IAt
        END IF
     END DO
  END DO


  if (debug) THEN
     WRITE(*,*) 'I found ',NbFrag, 'fragments with respectively '
     WRITE(*,*) (NbAtFrag(I),I=1,NbFrag), 'atoms'
     WRITE(*,*) "Fragments atoms are now listed as F in the following"
     DO I=1,NbFrag
        WRITE(*,*) Na
        WRITE(*,*) 'Fragment ', i
        DO J=1,Na
           AtName=Nom(Atome(J))
           IF (Fragment(J).EQ.I) AtName='F'
           WRITE(*,'(1X,A5,3(1X,F10.6))') AtName,X(J),Y(J),Z(J)
        END DO
        !     WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
     END DO

     DO I=1,NbFrag
        WRITE(*,*) 'Fragment ', i
        WRITE(*,'(A,20(1X,I3))') "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
        WRITE(*,'(A,20(1X,I3))') "FrozBlock:",(FrozBlock(I,j),j=0,FrozBlock(I,0))
     END DO
  END IF


  NFroz=0
  DO IFrag=1,NbFrag
     If (FrozBlock(IFrag,0).NE.0) NFroz=NFroz+1
  END DO

  IF (DEBUG) WRITE(*,'(1X,A,I4,A,I4,A)') "Found ",NFroz," fragmenst with frozen atoms, out of",NbFrag," fragments"

  ! We will now add the fragments containing non frozen atoms.
  ! I am not sure that there is a best strategy to add them...
  ! so we add them in the order they were created :-(
  ! First only block with frozen atoms are added
  izm0=-1
  IFrag=1
  FCaf=.FALSE.

  DO IFragFroz=1,NFroz
     DO WHILE (FrozBlock(IFrag,0).EQ.0)
        IFrag=IFrag+1
     END DO
     ! each atom has at least one frozen atom that will serve as the seed
     ! of this fragment.
     if (debug) WRITE(*,*) 'Adding fragment :', ifrag,'with',FrozBlock(IFrag,0),' frozen atoms'
     IF (FrozBlock(IFrag,0).EQ.1) THEN
        ! There is only one frozen atom, easy case ...
        Iat=FrozBlock(IFrag,1)
        if (debug) WRITE(*,*) 'Only frozen atom of frag',iat
        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)
           SELECT CASE(idx_zmat(n1))
           CASE (1)
              n2=ind_zmat(2,1)
              n3=ind_zmat(3,1)
           CASE (2)
              n2=ind_zmat(1,1)
              n3=ind_zmat(3,1)
           CASE (3:)
              n2=ind_zmat(idx_zmat(n1),2)
              n3=ind_zmat(idx_zmat(n1),3)
           END SELECT

           DO I3=1,Liaisons(n1,0)
              IAt=Liaisons(n1,I3)
              ! PFL 2007.Apr.16 ->
              ! We add ALL atoms linked to n1 to CaFaire
              ! But we delete their link to n1
              IF (.NOT.FCaf(Iat)) THEN
                 CaFaire(IdxCaFaire)=Iat
                 IdxCaFaire=IdxCaFaire+1
                 CaFaire(IdxCaFaire)=0
              END IF
              Call Annul(Liaisons,Iat,n1)
              Liaisons(iat,0)=Liaisons(Iat,0)-1
              ! <- PFL 2007.Apr.16
              IF (.NOT.DejaFait(Iat)) THEN
                 ! we add it to the Zmat
                 izm=izm+1
                 !              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)
                 idx_zmat(iat)=izm
                 !                  Call Annul(Liaisons,n1,iat)
                 n3=iat
                 DejaFait(Iat)=.TRUE.
              END IF
           END DO
              IaFaire=IaFaire+1

              if (debug) THEN
                 WRITE(*,*) 'ind_zmat 5'
                 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              ! DO WHILE CaFaire


        ELSE     ! FrozBlock(I,0)==1
           if (debug) WRITE(*,*) 'IFrag=',IFrag,'Frozblock(I,0)/=1',Frozblock(IFrag,0)
           ! We have many frozen atoms for one fragment. We will have to choose
           ! the first one, and then to decide when to swich...
           Iat=0
           I0=-1
           DO I=1,FrozBlock(IFrag,0)
              JAt=FrozBlock(IFrag,I)
              if (debug) WRITE(*,*) Jat, &
                   (LiaisonsBis(Jat,Itmp),Itmp=0,LiaisonsBis(Jat,0))
              IF (LiaisonsBis(Jat,0).GT.I0) THEN
                 I0=LiaisonsBis(Jat,0)
                 Iat=Jat
              END IF
           END DO
           ! Iat is the starting frozen atom
           IF (debug) WRITE(*,*) 'Choosing atom ',iat,'as a starting atom for frag',ifrag
           DejaFait(iat)=.TRUE.
           IdxCaFaire=2
           CaFaire(1)=iat
           CaFaire(2)=0
           IaFaire=1
           FCaf(Iat)=.TRUE.
           DO WHILE (CaFaire(IaFaire).NE.0)
              n1=CaFaire(IaFaire)
              if (debug) WRITE(*,*) 'Iafaire,n1,dejafait,idx_zmat', &
                   IaFaire,n1,    DejaFait(n1),idx_zmat(n1)
              if (debug) WRITE(*,*) 'Cafaire', &
                   (CaFaire(J),J=Iafaire,IdxCAfaire)
              SELECT CASE(idx_zmat(n1))
              CASE (1)
                 n2=ind_zmat(2,1)
                 n3=ind_zmat(3,1)
              CASE (2)
                 n2=ind_zmat(1,1)
                 n3=ind_zmat(3,1)
              CASE (3:)
                 n2=ind_zmat(idx_zmat(n1),2)
                 n3=ind_zmat(idx_zmat(n1),3)
              END SELECT

              if (debug) WRITE(*,*) "DBG n1,Liaisons(n1,0)",n1,Liaisons(n1,0)
              DO I3=1,Liaisons(n1,0)
                 IAt=Liaisons(n1,I3)                 
                 if (debug) WRITE(*,*) "DBG I3,Iat",I3,Iat
                 ! PFL 2007.Apr.16 ->
                 ! We add ALL atoms linked to n1 to CaFaire
                 ! But we delete their link to n1
!! PFL 2007.Aug.28 ->
! re-add the test on FCaf in order not to put the same atom more than once in 
! CaFaire array.
                 IF (.NOT.FCaf(Iat)) THEN
                    CaFaire(IdxCaFaire)=Iat
                    IdxCaFaire=IdxCaFaire+1
                    CaFaire(IdxCaFaire)=0
                    FCaf(Iat)=.TRUE.
                    if (debug) WRITE(*,*) 'Adding ',iat,' to CaFaire. I3=',I3,'IdxCafaire=',IdxCafaire
                    if (debug) WRITE(*,*) 'Cafaire',(CaFaire(J),J=Iafaire,IdxCAfaire)

                 END IF
!! <-- PFL 2007.Aug.28

                 Call Annul(Liaisons,Iat,n1)
                 Liaisons(iat,0)=Liaisons(Iat,0)-1
                 ! <- PFL 2007.Apr.16
                 IF (.NOT.DejaFait(Iat)) THEN
                    izm=izm+1
                    !                 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)
                    idx_zmat(iat)=izm
                    !                  Call Annul(Liaisons,n1,iat)

                    n3=iat
                    DejaFait(Iat)=.TRUE.
                 END IF
              END DO

              if (debug) WRITE(*,*) 'I3,n1,Liaisons(n1,:)', &
                   I3,n1,(Liaisons(n1,J),J=0,Liaisons(n1,0))


              if (debug.AND.(izm.GT.izm0)) THEN
                 WRITE(*,*) 'ind_zmat 6, 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 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
                 izm0=izm
              END IF

              IaFaire=IaFaire+1


           END Do              ! DO WHILE CaFaire

        END IF  ! FrozBlock(I,0)==1

        IFrag=IFrag+1
     END DO                    ! Loop on all fragments


     DO IFrag=1,NbFrag
        IF (FrozBlock(IFrag,0).EQ.0) THEN
           if (debug) WRITE(*,*) 'Adding fragment :', ifrag,'with no frozen atoms'
           ! We have no frozen atoms for this fragment. We will have to choose
           ! the first atom to put.
           ! For now, we do not care of the 'central' atom (ie the one with
           ! no CP atoms...)
           ! We just take the atom that is the closest to the already placed
           ! atoms !


           I0=0
           I1=0
           d=1e9
           DO J=1,izm
              Jat=ind_zmat(J,1)
              DO I=1,NbAtfrag(IFrag)
                 IAt=FragAt(IFrag,I)
                 IF (.NOT.DejaFait(Iat)) THEN
                    Call vecteur(Iat,Jat,x,y,z,vx1,vy1,vz1,norm1)
                    IF (norm1.LT.d) THEN
                       I1=Jat
                       I0=Iat
                       d=Norm1
                    END IF
                 END IF
              END DO
           END DO

           Iat=I0
           n1=I1

           IF (debug) WRITE(*,*) Iat,' starting atom for frag',ifrag


           SELECT CASE(idx_zmat(n1))
           CASE (1)
              n2=ind_zmat(2,1)
              n3=ind_zmat(3,1)
           CASE (2)
              n2=ind_zmat(1,1)
              n3=ind_zmat(3,1)
           CASE (3:)
              n2=ind_zmat(idx_zmat(n1),2)
              n3=ind_zmat(idx_zmat(n1),3)
           END SELECT

           izm=izm+1
           !        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)
           idx_zmat(iat)=izm


           DejaFait(iat)=.TRUE.
           IdxCaFaire=2
           CaFaire(1)=iat
           CaFaire(2)=0
           IaFaire=1
           FCaf(Iat)=.TRUE.
           DO WHILE (CaFaire(IaFaire).NE.0)
              n1=CaFaire(IaFaire)
              if (debug) WRITE(*,*) 'Iafaire,n1',IaFaire,n1, &
                   DejaFait(n1),idx_zmat(n1)
              if (debug) WRITE(*,*) 'Cafaire', &
                   (CaFaire(J),J=Iafaire,IdxCAfaire)
              SELECT CASE(idx_zmat(n1))
              CASE (1)
                 n2=ind_zmat(2,1)
                 n3=ind_zmat(3,1)
              CASE (2)
                 n2=ind_zmat(1,1)
                 n3=ind_zmat(3,1)
              CASE (3:)
                 n2=ind_zmat(idx_zmat(n1),2)
                 n3=ind_zmat(idx_zmat(n1),3)
              END SELECT


              DO I3=1,Liaisons(n1,0)
                 IAt=Liaisons(n1,I3)
                 ! PFL 2007.Apr.16 ->
                 ! We add ALL atoms linked to n1 to CaFaire
                 ! But we delete their link to n1
!! PFL 2007.Aug.28 ->
! re-add the test on FCaf in order not to put the same atom more than once in 
! CaFaire array.
                 IF (.NOT.FCaf(Iat)) THEN
                    CaFaire(IdxCaFaire)=Iat
                    IdxCaFaire=IdxCaFaire+1
                    CaFaire(IdxCaFaire)=0
                    FCaf(Iat)=.TRUE.
                 END IF
!! <-- PFL 2007.Aug.28

                 Call Annul(Liaisons,Iat,n1)
                 Liaisons(iat,0)=Liaisons(Iat,0)-1
                 if (debug) WRITE(*,*) 'Adding ',iat,' to CaFaire. I3=',I3
                 if (debug) WRITE(*,*) 'Cafaire',(CaFaire(J),J=Iafaire,IdxCAfaire)

                 ! <- PFL 2007.Apr.16
                 IF (.NOT.DejaFait(Liaisons(n1,i3))) THEN
                    IAt=Liaisons(n1,I3)
                    izm=izm+1
                    !                 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)
                    idx_zmat(iat)=izm

                    n3=iat
                    DejaFait(Iat)=.TRUE.
                 END IF

              if (debug) WRITE(*,*) 'I3,n1,Liaisons(n1,:)', &
                   I3,n1,(Liaisons(n1,J),J=0,Liaisons(n1,0))
              END DO
              IaFaire=IaFaire+1

              if (debug) THEN
                 WRITE(*,*) 'ind_zmat 7', 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 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              ! DO WHILE CaFaire
        END IF                 ! FrozBlock(IFrag,0).EQ.0

     END DO                    ! Loop on all fragments without frozen atoms

     if (debug) THEN
        WRITE(*,*) Na+Izm
        WRITE(*,*) 'Done... ', izm
        DO J=1,Na
           WRITE(*,'(1X,A5,3(1X,F10.6))') Nom(Atome(J)),X(J),Y(J),Z(J)
        END DO
        DO I=1,Izm
           J=ind_zmat(I,1)
           WRITE(*,'(1X,A5,3(1X,F10.6))') 'X ',X(J),Y(J),Z(J)
        END DO
        IF (izm.NE.Na) THEN
           WRITE(*,*) "Atoms not done:"
           DO I=1,Na
              IF (.NOT.DejaFait(I)) WRITE(*,'(I5)',ADVANCE='NO') I
           END DO
        END IF
     END IF


     ! 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,norm1,vx2,vy2,vz2,norm2, &
             vx4,vy4,vz4,norm4)
        CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
             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_Constr: ind_zmat'
        DO I=1,na
           WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
        END DO

        WRITE(*,*) 'DBG Cre_Zmat_Constr: 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 (debug) WRITE(*,*) "Deallocate (FrozDist,Fragment, NbAtFrag,FragAt)"
     DEALLOCATE(FrozDist,Fragment, NbAtFrag,FragAt)
     if (debug) WRITE(*,*) "Deallocate FrozFrag, IdxFragAt, FrozBlock"
     !  DEALLOCATE(FrozFrag, IdxFragAt, FrozBlock)
     DEALLOCATE(FrozFrag,FrozBlock)
     if (debug) WRITE(*,*) "Deallocate (DistFroz,Liaisons)"
     DEALLOCATE(DistFroz,Liaisons)
     if (debug) WRITE(*,*) "Deallocate (LiaisonsIni)"
     DEALLOCATE(LiaisonsIni)
     if (debug) WRITE(*,*) "Deallocate(CaFaire,DejaFait,FrozAt)"
     DEALLOCATE(CaFaire,DejaFait,FrozAt)
     if (debug) WRITE(*,*) "Deallocate (LiaisonsBis"
     DEALLOCATE(LiaisonsBis)


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

   END SUBROUTINE Calc_Zmat_const_frag


