      SUBROUTINE ZmatBuild(na,atome,ind_zmat,val_zmat,x,y,z,izm, &
           Liaisons, LIes_CP,LIes_CF,ListAt,DejaFait,debug)

!
!     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...
!
! 12.06.06
! Small modification: there was an inconsistency: izm was increased
! after the atom was added to ind_zmat for the three first atoms,
! but before for all others. Now, it starts at zero and it
! is increased before the atom is added to ind_zmat
!     IMPLICIT NONE

      use Path_module, only : Nat, NMaxL, max_Z, KINT, KREAL


      INTEGER(KINT) :: Izm
      integer(KINT) :: na,atome(NAt),at,ind_zmat(NAt,5),long
      real(KREAL) ::  x(NAt),y(NAt),z(NAt),fact
      real(KREAL) ::  val_zmat(NAt,3)
!     ListAt contains the indices of frozen atoms
      LOGICAL ListAT(NAt),FIzm1,FFirst
      INTEGER(KINT) :: Natreat

      real(KREAL) ::  vx,vy,vz,dist
      INTEGER(KINT) :: LIAISONS(NAt,0:NMaxL),Nbli,Nblj
      INTEGER(KINT) :: LiaisonsBis(NAt,0:NMaxL)
      INTEGER(KINT) :: CAFaire(NAt)
      INTEGER(KINT) :: LieS_CP(NAt,0:NMaxL),LieS_CF(NAt,0:NMaxL)
      INTEGER(KINT) :: AtCP0(NAt),NbAt0,NbAt0r
      Integer(KINT) :: AtTypCycl(Max_Z), NbAt3(NAt,2),CyCleAt(NAt)
      Integer(KINT) :: NbCycle
      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
      Logical AtTerm,Debug,Done
      Logical DejaFait(NAt)
      LOGICAL F1213, F1223,F1323,FTmp

!     As we may have to treat only partial molecules, it may happen
!     that there are no central atoms... so we first check this by looking
!     for the least number of CP atoms
      FFirst=.TRUE.
      NMinCP=na
      NMaxCF=-1
      NAtreat=0
      Izm0=Izm+1

      DO I=1,na
         if (DEBUG) WRITE(*,*) "DBG ZmatBuild i,listat",i,listAt(i)
         IF (ListAt(I).AND.(.NOT.(DejaFait(I)))) THEN
            IF (Lies_CP(I,0).lt.NMinCP)      NMinCP=Lies_CP(I,0)
            IF (Lies_CF(I,0).gt.NMaxCF)      NMaxCF=Lies_CF(I,0)
            NATreat=NATreat+1
         END IF
      END DO
      if (debug) WRITE(*,*) 'NMinCP=',NMinCP
      if (debug) WRITE(*,*) 'NMaxCF=',NMaxCF
      if (debug) WRITE(*,*) 'NATreat=',NAtreat


      IF (Debug) THEN
         WRITE(*,*) '------------------ Zmat Build -------------------'
         WRITE(*,*) 'LIAISONS'
         DO I=1,na
            WRITE(*,1003) i,(LIAISONS(I,J),J=0,NMaxL)
         END DO

         WRITE(*,*) 'LIes_CF'
         DO I=1,na
            WRITE(*,1003) i,(LIeS_CF(I,J),J=0,NMaxL)
         END DO

         WRITE(*,*) 'LIes_CP'
         DO I=1,na
            WRITE(*,1003) i,(LIeS_CP(I,J),J=0,NMaxL)
         END DO
         WRITE(*,*) '-- Zmat Build : only _To treat_ atoms------------'
         WRITE(*,*) 'LIAISONS'
         DO I=1,na
            IF (ListAt(I).AND.(.NOT.DejaFait(I))) &
                 WRITE(*,1003) i,(LIAISONS(I,J),J=0,NMaxL)
         END DO

         WRITE(*,*) 'LIes_CF'
         DO I=1,na
            IF (ListAt(I).AND.(.NOT.DejaFait(I))) &
                WRITE(*,1003) i,(LIeS_CF(I,J),J=0,NMaxL)
         END DO

         WRITE(*,*) 'LIes_CP'
         DO I=1,na
            IF (ListAt(I).AND.(.NOT.DejaFait(I))) &
                WRITE(*,1003) i,(LIeS_CP(I,J),J=0,NMaxL)
         END DO

      END IF


      if (NAtreat.EQ.0) THEN
         WRITE(*,*) "Foutage de gueule : NAtTreat=0"
         RETURN
      END IF

      if (debug) THEN
         WRITE(*,*) 'DBG ZmatBuil, izm=',izm
         DO I=1,na
            WRITE(*,'(1X,I5,1X,L4,1X,L4,12(1X,I4))') i,ListAt(I) &
                 ,DejaFait(I), &
                 (Liaisons(I,J),J=0,Liaisons(I,0))
         END DO
      END IF

!!!   atraiter NMaxCF=0 :que des atomes separes...
!     mais faut-il reellement faire un cas a part ?

!     On compte le nb d'atomes sans atomes CP (centripetes)
      NbAt0=0
      DO I=1,na
         IF (ListAt(I).AND.(.NOT.DejaFait(I)).AND. &
              (Lies_CP(I,0).eq.NMinCP)) THEN
            NbAt0=NbAt0+1
            AtCP0(NbAt0)=I
            IF (DEBUG) WRITE(*,*) "Center atom",NbAt0,I,Atome(I)
         END IF
      END DO
      AtCP0(NbAt0+1)=0
      NbAt0r=NbAt0

      IF (Debug)  WRITE(*,*) 'DBG ZmatBuld - NCP0,AtCP0 ',NbAt0, &
           (AtCP0(i),i=1,NbAt0)

!     On va les traiter tous en construisant les molecules en partant d'eux
!     Le premier atome sans CP est different des autres car il fixe
!     l'origine

!First atom
      IF (Izm==0) THEN
         FFirst=.FALSE.
         if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=1'
!     IZm=1
!     Boucle pour trouver celui qui a le plus d'atomes CF
         IdAt0=0
         VCF=-1
         DO I=1,NbAt0
!     WRITE(*,*) 'DBG ZmatB',I,AtCP0(I)
            if (AtCP0(I).NE.0) THEN
               ICat=AtCP0(I)
               if (Lies_CF(ICAt,0).gt.VCF) THEN
                  IdAt0=ICat
                  IdxAt0=I
                  VCF=Lies_CF(ICAt,0)
               END IF
            END IF
         END DO

!     IF (debug) WRITE(*,*) 'DBG ZmatBuil - VCF, IdxAt0,IdAt0',
!     &        VCF, IdxAt0,IdAt0
!     all atoms with NMinCP CP links do not have CF links... not a good choice to start
!     building the molecule(s) around them... we increase NMinCP
         AtCP0(IdxAt0)=0
         NbAt0r=NbAt0r-1

!     WRITE(IOOUT,*) 'Atome 1', IdAt0, Nom(Atome(IdAt0))
         Izm1=IdAt0
         Izm=Izm+1

         ind_zmat(izm,1)=Izm1
         ind_zmat(izm,2)=0
         ind_zmat(izm,3)=0
         ind_zmat(izm,4)=0
         ind_zmat(izm,5)=0
         val_zmat(izm,1)=0.0
         val_zmat(izm,2)=0.0
         val_zmat(izm,3)=0.0
         DejaFait(Izm1)=.True.

      END IF                    ! end of izm==1 test

!     On construit la premiere couche qui est speciale car elle fixe les
!     axes.
!     Plusieurs cas sont possibles suivant le nb d'atomes CF
      IdAt0=ind_zmat(1,1)
      Izm1=IdAt0

!Second atom
      IF ((izm==1).AND.(NAtreat.ge.2))  THEN
         Done=.FALSE.
         if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=2'
!     Two cases: 1) We are constructing the first fragment: FIzm1=.T.
!     2) We are constructing the second fragment: FIzm1=.F.
!     in this case, we select one CP0 atom to start with
         IF (FFirst) THEN
            FFirst=.FALSE.
!     This is the first atom to be dealt with
!     We look for a CP0
            If (NbAt0r.ge.1) THEN
               IdAt0=0
               VCF=-1
               WRITE(*,*) 'DBG  ZBuild Izm=2:AtCP0',(AtCP0(I),I=1,NbAt0)
               DO I=1,NbAt0
                  if (AtCP0(I).NE.0) THEN
                     ICat=AtCP0(I)
                     if (Lies_CF(I,0).gt.VCF) THEN
                        Izm2=ICat
                        IdxAt0=I
                        VCF=Lies_CF(ICAt,0)
                     END IF
                  END IF
               END DO
               WRITE(*,*) 'DBG ZBuild Izm=2:izm2,IdxAt0,VCF,AtCP0(1)', &
                    izm2,IdxAt0,VCF,AtCP0(1)
               NbAt0r=NbAt0r-1
               Done=.TRUE.
            END IF

!     If we do not have a CP0 available the other tests are identical
!     for cases 1 and 2...
         END IF
         IF (.NOT.DONE) THEN
            IF (Lies_CF(IdAt0,0).ge.1) THEN
               IZm2=Lies_CF(IdAt0,1)
               WRITE(*,*) 'DBG ZBuild Lies_CF(IdAt0,0).ge.1', &
                    Lies_CF(IdAt0,0),izm2
            ELSE
!     First atom is not linked to anything... we look for the second CP0 atom
!     if we have one..
               If (NbAt0r.ge.1) THEN
                  IdAt0=0
                  VCF=-1
                  IF (DEBUG) WRITE(*,*) 'DBG  ZBuild Izm=2:AtCP0', &
                       (AtCP0(I),I=1,NbAt0)
                  DO I=1,NbAt0
                     if (AtCP0(I).NE.0) THEN
                        ICat=AtCP0(I)
                        if (Lies_CF(I,0).gt.VCF) THEN
                           Izm2=ICat
                           IdxAt0=I
                           VCF=Lies_CF(ICAt,0)
                        END IF
                     END IF
                  END DO
                  IF (DEBUG) WRITE(*,*) 'DBG ZBuild Izm=2:izm2,IdxAt0,VCF,AtCP0(1)', &
                       izm2,IdxAt0,VCF,AtCP0(1)
                  NbAt0r=NbAt0r-1
               ELSE
!     we do not have another CP0 atom... but we know that there
!     is at least two atoms... we look for the closest atom to Izm1
                  Izm2=0
                  Dist=1.D99
                  DO I=1,Na
                     IF (ListAt(I).AND.(.NOT.Dejafait(I))) THEN
                        CALL vecteur(IZm1,I,x,y,z,vx1,vy1,vz1,norm1)
                        if (Norm1.lt.Dist) IZm2=I
                     END IF
                  END DO
               END IF
            END IF
         END IF
         CALL vecteur(IZm1,IZm2,x,y,z,vx1,vy1,vz1,norm1)

         Izm=Izm+1

         ind_zmat(izm,1)=izm2
         ind_zmat(izm,2)=izm1
         ind_zmat(izm,3)=0
         ind_zmat(izm,4)=0
         ind_zmat(izm,5)=0
         val_zmat(izm,1)=norm1
         val_zmat(izm,2)=0.0
         val_zmat(izm,3)=0.0
         DejaFait(Izm2)=.TRUE.

      END IF                    ! end of izm==2 test

      Izm1=ind_zmat(1,1)
      Izm2=ind_zmat(2,1)

! Third atom

      IF ((Izm==2).AND.(NAtreat.ge.3)) THEN
         Done=.FALSE.
         if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=3',FFirst
!     Two cases: 1) We are constructing the first fragment: FIzm1=.T.
!     2) We are constructing the second fragment: FIzm1=.F.
!     in this case, we select one CP0 atom to start with
         IF (FFirst) THEN
            FFirst=.FALSE.
!     This is the first atom to be dealt with
!     We look for a CP0
            If (NbAt0r.ge.1) THEN
               IdAt0=0
               VCF=-1
               DO I=1,NbAt0
                  if (AtCP0(I).NE.0) THEN
                     ICat=AtCP0(I)
                     if (Lies_CF(I,0).gt.VCF) THEN
                        Izm3=ICat
                        IdxAt0=I
                        VCF=Lies_CF(ICAt,0)
                     END IF
                  END IF
               END DO
               NbAt0r=NbAt0r-1
               Done=.TRUE.
            END IF
!     If we do not have a CP0 available the other tests are identical
!     for cases 1 and 2...
         END IF
         IF (.NOT.Done) THEN
            WRITE(*,*) 'DBG ZmatBuild: Done false for izm=3'
            IF (Lies_CF(izm1,0).ge.2) THEN
!     easiest case: the first atom is linked to at least 2 atoms.
               WRITE(*,*) 'DBG ZmtB-izm3:Lies_CF(izm1,0)>=2 '
               I=1
               DO WHILE ((.NOT.ListAt(Lies_CF(izm1,I))) &
                    .OR.(DejaFait(Lies_CF(izm1,I))))
                  I=I+1
               END DO
               izm3=Lies_CF(Izm1,I)
               Done=.TRUE.
            ELSEIF (Lies_CF(Izm2,0).ge.1) THEN
!     a bit more complex: first atom is linked to one atom (Izm2)
!     but luckily 2nd atom is linked to at least one atom
               I=1
               DO WHILE ((.NOT.ListAt(Lies_CF(izm2,I))) &
                    .OR.(DejaFait(Lies_CF(izm2,I))).AND.(I.le.Lies_CF(Izm2,0)))
                  I=I+1
               END DO
               IF (I.LE.Lies_CF(Izm2,0)) THEN
                  Izm3=Lies_CF(izm2,I)
               WRITE(*,*) 'DBG ZmtB-izm3:Lies_CF(izm2,0)>=1 and ok  '
!     We exchange Izm1 and Izm2 because the logical order is Izm3 linked
!     to Izm2  and not to Izm1 as is the default later
                  IzmTmp=Izm2
                  Izm2=Izm1
                  Izm1=IzmTmp
                  Done=.TRUE.
               END IF
            END IF
         IF (.NOT.Done) THEN
!     None of the first two atoms has links to a third atom...
!     we take it from the CP0 if possible...
               If (NbAt0r.ge.1) THEN
                  WRITE(*,*) 'DBG ZmatB izm=3',NbAt0r
                  IdAt0=0
                  VCF=-1
                  DO I=1,NbAt0
                     if (AtCP0(I).NE.0) THEN
                        ICat=AtCP0(I)
                        if (Lies_CF(I,0).gt.VCF) THEN
                           Izm3=ICat
                           IdxAt0=I
                           VCF=Lies_CF(ICAt,0)
                        END IF
                     END IF
                  END DO
                  NbAt0r=NbAt0r-1
                  WRITE(*,*) 'DBG ZmatB izm=3',NbAt0r,Izm3
               ELSE
! Or the  atom closest to 1
                  Izm3=0
                  Dist=1.D99
                  DO I=1,Na
                     IF (ListAt(I).AND.(.NOT.Dejafait(I))) THEN
                        CALL vecteur(IZm1,I,x,y,z,vx1,vy1,vz1,norm1)
                        if (Norm1.lt.Dist) IZm3=I
                     END IF
                  END DO
               END IF
            END IF
         END IF

         CALL vecteur(IZm1,IZm2,x,y,z,vx1,vy1,vz1,norm1)
         CALL vecteur(izm1,izm3,x,y,z,vx2,vy2,vz2,norm2)
         val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)

!     write(*,11) n1,n2,norm1,n3,val

         Izm=Izm+1

         ind_zmat(izm,1)=izm3
         ind_zmat(izm,2)=izm1
         ind_zmat(izm,3)=izm2
         ind_zmat(izm,4)=0
         ind_zmat(izm,5)=0
         val_zmat(izm,1)=norm2
         val_zmat(izm,2)=val
         val_zmat(izm,3)=0.0
         DejaFait(Izm3)=.TRUE.

      END IF                    ! end of test izm=3

!     We add all atoms CF atoms linked to the already present atoms in the zmat to
!     the "to do" list:
!     Les atomes CF lies a IdAt0 sont mis en attente pour
!     etre traites

!     For 'historical' reasons, the first atom links have to be dealt with
!     in a special way... now !

      if (Izm.ge.NaTreat) Return


      WRITE(*,*) 'DBG ZmatBuild Izm0',Izm0

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

      IF (FFIrst) GOTO 9753

         I=Ind_zmat(Izm0,1)
         DO IdAtl=1,Lies_CF(I,0)
            Izm4= Lies_CF(I,IdAtl)
!      WRITE(IOOUT,*) 'Atome ',Izm+1,Izm4,Nom(Atome(Izm4)),Izm1,
!     $Izm2,Izm3
            IF (ListAt(Izm4).AND.(.NOT.DejaFait(Izm4))) THEN
               Izm=Izm+1
               Call AddLigne(Izm,Izm4,Izm1,Izm2,Izm3,ind_zmat,val_zmat, &
                 x,y,z)
               DejaFait(Izm4)=.TRUE.
            END IF
         END DO
!      END IF

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

      if (Izm.ge.NaTreat) Return

      IndFin=1
      DO I=Izm0,Izm
         Iat=ind_zmat(I,1)
!         Do iaf=1,Lies_CF(Iat,0)
!!     if (ListAt(Lies_CF(IAt,iaf)).AND.
!!     &              (.NOT.DejaFait(Lies_CF(IAt,iaf))))  THEN
!            if (ListAt(Lies_CF(IAt,iaf))) THEN
!               CAfaire(IndFin)=Lies_CF(IAt,Iaf)
               CaFaire(IndFin)=Iat
               IndFin=IndFin+1
!            END IF
!         END DO
      END DO
      CAfaire(IndFin)=0

      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire

      if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=4'

!     on a cree la premiere couche autour du premier centre
!     reste a finir les autres couches
      IndAFaire=1
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
      Do WHILE (Cafaire(IndAFaire).ne.0)
         IAt0=Cafaire(IndAfaire)
         if (debug) WRITE(*,*) "IndAFaire,IAt0,Lies_CF(IAt0,0)", &
              IndAFaire,IAt0,Lies_CF(IAt0,0)
         if (Lies_CF(IAt0,0).ge.1) THEN
!     IAt0 est l'atome pour lequel on construit la couche suivante
!     le premier atome lie est particulier car il definit l'orientation
!     de ce fragment
            IAti=Lies_CF(IAt0,1)
      WRITE(IOOUT,*) 'IAti:',IAti
            IAtd=IAt0
      WRITE(IOOUT,*) 'IAtd:',IAtd
            IAtv=Lies_CP(IAt0,1)
      WRITE(IOOUT,*) 'IAtv:',IAtv
            IIAtdh=1
            IAtdh=Lies_CF(IAtv,IIatdh)
            DO WHILE (Iatdh.eq.IAt0)
               IIatdh=IIatdh+1
               IAtdh=Lies_CF(IAtv,IIatdh)
            END DO
            IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
      WRITE(IOOUT,*) 'IAtdh:',IAtdh
            IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
               Izm=Izm+1
               if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
                    izm,IAti,IAtd,IAtv,IAtdh
               Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
                    ind_zmat,val_zmat ,x,y,z)
               DejaFait(IAti)=.TRUE.
            END IF

            CAfaire(IndFin)=Lies_CF(IAt0,1)
            IndFin=IndFin+1
!     Le traitement des autres est plus facile
            IAtdh=Lies_CF(IAt0,1)
            DO IAta=2,Lies_CF(IAt0,0)
               IAtI=Lies_CF(IAt0,IAta)
               CAfaire(IndFin)=IAtI
               IndFin=IndFin+1

               IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
                  Izm=Izm+1
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh,ind_zmat &
                       ,val_zmat,x,y,z)
                  DejaFait(IAti)=.TRUE.
               END IF
            END DO
            CAfaire(IndFin)=0
         END IF
         IndAFaire=IndAFaire+1
         if (debug) WRITE(IOOUT,*) "IndAfaire,CaFaire:",IndAfaire, &
            CaFaire
      END DO

      IndFin=1
      DO I=1,Izm0
         Iat=ind_zmat(I,1)
!         Do iaf=1,Lies_CF(Iat,0)
         if (ListAt(IAt).AND. &
                    (.NOT.DejaFait(IAt)))  THEN
!            if (ListAt(Lies_CF(IAt,iaf))) THEN
!               CAfaire(IndFin)=Lies_CF(IAt,Iaf)
               CaFaire(IndFin)=Iat
               IndFin=IndFin+1
            END IF
!         END DO
      END DO
      CAfaire(IndFin)=0

      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire

      if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=4'

!     on a cree la premiere couche autour du premier centre
!     reste a finir les autres couches
      IndAFaire=1
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
      Do WHILE (Cafaire(IndAFaire).ne.0)
         IAt0=Cafaire(IndAfaire)
         if (debug) WRITE(*,*) "IndAFaire,IAt0,Lies_CF(IAt0,0)", &
              IndAFaire,IAt0,Lies_CF(IAt0,0)
         if (Lies_CF(IAt0,0).ge.1) THEN
!     IAt0 est l'atome pour lequel on construit la couche suivante
!     le premier atome lie est particulier car il definit l'orientation
!     de ce fragment
            IAti=Lies_CF(IAt0,1)
      WRITE(IOOUT,*) 'IAti:',IAti
            IAtd=IAt0
      WRITE(IOOUT,*) 'IAtd:',IAtd
            IAtv=Lies_CP(IAt0,1)
      WRITE(IOOUT,*) 'IAtv:',IAtv
            IIAtdh=1
            IAtdh=Lies_CF(IAtv,IIatdh)
            DO WHILE (Iatdh.eq.IAt0)
               IIatdh=IIatdh+1
               IAtdh=Lies_CF(IAtv,IIatdh)
            END DO
            IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
      WRITE(IOOUT,*) 'IAtdh:',IAtdh
            IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
               Izm=Izm+1
               if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
                    izm,IAti,IAtd,IAtv,IAtdh
               Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
                    ind_zmat,val_zmat ,x,y,z)
               DejaFait(IAti)=.TRUE.
            END IF

            CAfaire(IndFin)=Lies_CF(IAt0,1)
            IndFin=IndFin+1
!     Le traitement des autres est plus facile
            IAtdh=Lies_CF(IAt0,1)
            DO IAta=2,Lies_CF(IAt0,0)
               IAtI=Lies_CF(IAt0,IAta)
               CAfaire(IndFin)=IAtI
               IndFin=IndFin+1

               IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
                  Izm=Izm+1
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh,ind_zmat &
                       ,val_zmat,x,y,z)
                  DejaFait(IAti)=.TRUE.
               END IF
            END DO
            CAfaire(IndFin)=0
         END IF
         IndAFaire=IndAFaire+1
         if (debug) WRITE(IOOUT,*) "IndAfaire,CaFaire:",IndAfaire, &
            CaFaire
      END DO


!     On a fini de creer la molecule autour du premier atome 'centre'
!     Pour les autres c'est presque id... sauf que les axes sont deja fixes
!     On pourrait imaginer de mettre des atomes fictifs pour redefinir des axes
!     locaux... dans une autre version
!     V2.0
 9753 FFirst=.FALSE.
      if (debug) THEN
         WRITE(*,*) 'NbAt0r',NbAt0r,AtCP0
      END IF
      DO I=1,NbAt0r
!     Boucle pour trouver celui qui a le plus d'atomes CF
         IdAt0=0
         VCF=-1
         DO ICP0=1,NbAt0
            if (AtCP0(ICP0).NE.0) THEN
               ICat=AtCP0(ICP0)
               if (Lies_CF(ICAt,0).gt.VCF) THEN
                  IdAt0=ICat
                  IdxAt0=ICP0
                  VCF=Lies_CF(ICAt,0)
               END IF
            END IF
         END DO
         AtCP0(IdxAt0)=0

         IF (debug) WRITE(*,*) 'Adding fragment centered on ',IdAt0, &
              LIAISONS(IdAt0,0)
         IF (debug) WRITE(*,*) 'Center linked to ',Lies_CF(IdAt0,0), &
              ' atoms'
!     if (LIAISONS(IdAt0,0).EQ.0) goto 12345

         if (debug) THEN
            WRITE(IOOUT,*) 'LIAISONS tralala',nbcycle
            DO IAt=1,na
               WRITE(IOOUT,1003) iat,(LIAISONS(Iat,Jat),Jat=0,NMaxL)
            END DO
         END IF
 1003    FORMAT(1X,I5,' - ',12(I5))
!     Boucle pour voir quel est l'atome du fragment precedent qui est le plus
!     proche de celui ci
!     Rq: si ce sont des cycles lies, il pourrait etre malin de rechercher
!     a quoi il etait lie.
         norm2=1.e6
         Idx=0
         DO ICAt=1,Izm
            CALL vecteur(Idat0,ind_zmat(icat,1),x,y,z,vx1,vy1,vz1 &
                 ,norm1)
            if (norm2.ge.norm1) THEN
               norm2=norm1
               Idx=Icat
            END IF
         END DO
         IF (debug) WRITE(*,*) 'Link between fragments ',IdAt0, &
              ind_zmat(Idx,1), ' -- Idx=',Idx

! Despite the fact that this atom has officialy no CP atom
! we add one into its list; the one just found
         Lies_CP(IdAt0,0)=Lies_CP(IdAt0,0)+1
         Lies_CP(Idat0,Lies_CP(IdAt0,0))= ind_zmat(Idx,1)
         Lies_CP(Idat0,Lies_CP(IdAt0,0)+1)=0


!     on a l'indice... on va rajouter cet atome a la liste !
         izm=izm+1
         n1=ind_zmat(Idx,1)
!     Petit probleme si Idx<=2
         IF (Idx.EQ.1) THEN
!     Pb non negligeable ...
            IF (izm.ge.2) THEN
               IAtv=ind_zmat(2,1)
               IAtdh=ind_zmat(3,1)
            ELSE
               WRITE(*,*) '2 frag, le 1er a moins de 2 atomes..'
               WRITE(*,*) "J'ai merde... "
               STOP
            END IF
         ELSEIF (Idx.EQ.2) THEN
            IAtv=1
            IF (izm.ge.2) THEN
               IAtdh=ind_zmat(3,1)
            ELSE
               WRITE(*,*) '2 frag, le 1er a moins de 2 atomes...'
               WRITE(*,*) "J'ai merde... "
               STOP
            END IF
         ELSE
            IAtv=ind_zmat(Idx,2)
            IAtdh=ind_zmat(Idx,3)
         END IF

!     WRITE(*,*) "Outchilou IAtv, IAtdh", IAtv, IAtdh,n1
         CALL vecteur(n1,IdAt0,x,y,z,vx1,vy1,vz1,norm1)
         CALL vecteur(n1,IAtv,x,y,z,vx2,vy2,vz2,norm2)
         val=angle(vx1,vy1,vz1,norm1, &
              vx2,vy2,vz2,norm2)
         if (debug) WRITE(*,*) 'Angle val:1',IdAt0,n1,IAtv,val
         if ((abs(val).LE.10.).OR.(abs(pi-val).le.10.)) THEN
            IAtv=IAtdh
!     Ceci pose pb si Idx<=3... a traiter plus tard
            IF (Idx.ge.3) THEN
               IAtdh=ind_zmat(Idx,4)
            ELSE
               if (izm.ge.4) THEN
                  IAtdh=4
               ELSE
                  WRITE(*,*) 'Pb Angle Val+idx=3+izm<=4'
                  STOP
               END IF
            END IF
         END IF
         Call AddLigne(Izm,IdAt0,ind_zmat(Idx,1),IAtv,IAtdh, &
              ind_zmat,val_zmat,x,y,z)
         DejaFait(IdAt0)=.TRUE.

         IndFin=1
         IAtd=IdAt0
         n1=IAtdh
         IAtdh=IAtv
         IAtv=ind_zmat(Idx,1)
!     On ajoute les atomes lies a cet atome dans la liste a faire
         Do iaf=1,Lies_CF(IdAt0,0)
            IatI=Lies_CF(IdAt0,Iaf)
!     We check that this atom is not the one that is the closest to
!     the center atom
            if (IAtv.Ne.IAtI) THEN
               IF (debug) WRITE(*,*) 'Adding atom',IAtI, &
                    'to CAFaire in pos',iaf
               CAfaire(IndFin)=IAtI
               IndFin=IndFin+1
!     On les ajoute aussi dans la zmat
               If (ListAt(IAtI).AND.(.NOT.DejaFait(IatI))) THEN
                  izm=izm+1
       WRITE(*,*) "Outchili IAtv, IAtdh,IAtI", IAtv, IAtdh,IAtI
                  CALL vecteur(IAtI,IAtD,x,y,z,vx1,vy1,vz1,norm1)
                  CALL vecteur(IAtI,IAtv,x,y,z,vx2,vy2,vz2,norm2)
                  val=angle(vx1,vy1,vz1,norm1, vx2,vy2,vz2,norm2)
                  if (debug) &
                       WRITE(*,*) 'Angle val:2',IAtI,IAtD,IAtv,val
                  if ((abs(val).LE.10.).OR. &
                       (abs(180.-abs(val)).le.10.)) THEN
                     Call AddLigne(Izm,IAtI,IAtv,IAtdh,n1, &
                          ind_zmat,val_zmat,x,y,z)
                     DejaFait(Iati)=.TRUE.
                  ELSE
                     Call AddLigne(Izm,IAtI,IAtD,IAtv,IAtdh, &
                          ind_zmat,val_zmat,x,y,z)
                     DejaFait(Iati)=.TRUE.
                  END IF
               END IF
            END IF
         END DO
         CAfaire(IndFin)=0

!     On traite le reste de ce fragment !!
         IndAFaire=1
         WRITE(IOOUT,*) CaFaire
         Do WHILE (Cafaire(IndAFaire).ne.0)
            IAt0=Cafaire(IndAfaire)
            if (Lies_CF(IAt0,0).ge.1) THEN
!     IAt0 est l'atome pour lequel on construit la couche suivante
!     le premier atome lie est particulier car il definit l'orientation
!     de ce fragment
               Itmp=1
               IAti=Lies_CF(IAt0,Itmp)
               DO WHILE (DejaFait(IAti).AND. &
                    (Itmp.Le.Lies_CF(IAt0,0)))
                  Itmp=Itmp+1
                  IAti=Lies_CF(IAt0,ITmp)
               END DO
               If (ListAt(IAtI).AND.(.NOT.DejaFait(IatI))) THEN
                  IF (debug) WRITE(IOOUT,*) 'IAti:',IAti
                  IAtd=IAt0
                  IF (debug)     WRITE(IOOUT,*) 'IAtd:',IAtd
                  IAtv=Lies_CP(IAt0,1)
                  IF (debug)     WRITE(IOOUT,*) 'IAtv:',IAtv
                  IIAtdh=1
                  IAtdh=Lies_CF(IAtv,IIatdh)
                  DO WHILE (Iatdh.eq.IAt0)
                     IIatdh=IIatdh+1
                     IAtdh=Lies_CF(IAtv,IIatdh)
                  END DO
                  IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
                  IF (debug) WRITE(IOOUT,*) 'IAtdh:',IAtdh
                  Izm=Izm+1
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
                       ind_zmat,val_zmat,x,y,z)
                     DejaFait(Iati)=.TRUE.
!     Le traitement des autres est plus facile
                  IAtdh=Lies_CF(IAt0,ITmp)
                  DO IAta=ITmp+1,Lies_CF(IAt0,0)
                     IAtI=Lies_CF(IAt0,IAta)
                     If (ListAt(IAtI).AND.(.NOT.DejaFait(IatI))) &
                          THEN
                        Izm=Izm+1
                        Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
                             ind_zmat,val_zmat,x,y,z)
                        DejaFait(Iati)=.TRUE.
                     END IF
                  END DO
               END IF

!     On ajoute les atomes lies a cet atome dans la liste a faire
               Do iaf=1,Lies_CF(IAt0,0)
                  CAfaire(IndFin)=Lies_CF(IAt0,Iaf)
                  IndFin=IndFin+1
               END DO
               CAfaire(IndFin)=0
            END IF
            IndAFaire=IndAFaire+1
         END DO
12345    CONTINUE
      END DO

      WRITE(*,*) 'TOTOTOTOTOTOTOTOT'

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

      IndFin=1
      DO I=1,Izm0
         Iat=ind_zmat(I,1)
         Do iaf=1,Lies_CF(Iat,0)
            IIat=Lies_CF(Iat,iaf)
            if (ListAt(iIAt).AND. &
                 (.NOT.DejaFait(iIAt)))  THEN
               IAti=IIat
               WRITE(IOOUT,*) 'IAti:',IAti
               IAtd=IAt0
               WRITE(IOOUT,*) 'IAtd:',IAtd
               IAtv=Lies_CP(IAt0,1)
               WRITE(IOOUT,*) 'IAtv:',IAtv
               IIAtdh=1
               IAtdh=Lies_CF(IAtv,IIatdh)
               DO WHILE (Iatdh.eq.IAt0)
                  IIatdh=IIatdh+1
                  IAtdh=Lies_CF(IAtv,IIatdh)
               END DO
               IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
               WRITE(IOOUT,*) 'IAtdh:',IAtdh
               IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
                  Izm=Izm+1
                  if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
                       izm,IAti,IAtd,IAtv,IAtdh
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
                       ind_zmat,val_zmat ,x,y,z)
                  DejaFait(IAti)=.TRUE.
                  CaFaire(IndFin)=iIat
                  IndFin=IndFin+1
               END IF
            END IF
         END DO
      END DO
      CAfaire(IndFin)=0

      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire

      if (debug) WRITE(*,*) 'DBG ZmatBuil - Toto'
      if (debug) THEN
         WRITE(*,*) 'ind_zmat'
         DO I=1,na
            WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
         END DO
      END IF

!     on a cree la premiere couche autour du premier centre
!     reste a finir les autres couches
      IndAFaire=1
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
      Do WHILE (Cafaire(IndAFaire).ne.0)
         IAt0=Cafaire(IndAfaire)
         if (debug) WRITE(*,*) "IndAFaire,IAt0,Lies_CF(IAt0,0)", &
              IndAFaire,IAt0,Lies_CF(IAt0,0)
         if (Lies_CF(IAt0,0).ge.1) THEN
!     IAt0 est l'atome pour lequel on construit la couche suivante
!     le premier atome lie est particulier car il definit l'orientation
!     de ce fragment
            IAti=Lies_CF(IAt0,1)
      WRITE(IOOUT,*) 'IAti:',IAti
            IAtd=IAt0
      WRITE(IOOUT,*) 'IAtd:',IAtd
            IAtv=Lies_CP(IAt0,1)
      WRITE(IOOUT,*) 'IAtv:',IAtv
            IIAtdh=1
            IAtdh=Lies_CF(IAtv,IIatdh)
            DO WHILE (Iatdh.eq.IAt0)
               IIatdh=IIatdh+1
               IAtdh=Lies_CF(IAtv,IIatdh)
            END DO
            IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
      WRITE(IOOUT,*) 'IAtdh:',IAtdh
            IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
               Izm=Izm+1
               if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
                    izm,IAti,IAtd,IAtv,IAtdh
               Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
                    ind_zmat,val_zmat ,x,y,z)
               DejaFait(IAti)=.TRUE.
            END IF

            CAfaire(IndFin)=Lies_CF(IAt0,1)
            IndFin=IndFin+1
!     Le traitement des autres est plus facile
            IAtdh=Lies_CF(IAt0,1)
            DO IAta=2,Lies_CF(IAt0,0)
               IAtI=Lies_CF(IAt0,IAta)
               CAfaire(IndFin)=IAtI
               IndFin=IndFin+1

               IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
                  IF (debug) &
         WRITE(*,*) 'DBG ZmatBuil Toto adding',Iati,Izm
                  Izm=Izm+1
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh,ind_zmat &
                       ,val_zmat,x,y,z)
                  DejaFait(IAti)=.TRUE.
               END IF
            END DO
            CAfaire(IndFin)=0
         END IF
         IndAFaire=IndAFaire+1
         if (debug) WRITE(IOOUT,*) "Toto IndAfaire,CaFaire:",IndAfaire, &
            CaFaire
      END DO


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

      END
