SUBROUTINE Calc_mixed_frag(na,atome,ind_zmat,val_zmat,x,y,z, &
     r_cov,fact,frozen,cart,ncart)

  !     
  !     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(*),Cart(*),NFroz,NCart
  LOGICAL, ALLOCATABLE :: FCart(:) ! Na
  INTEGER(KINT), ALLOCATABLE :: AtCart(:) !Na
  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,IBeg
  INTEGER(KINT) :: I0, Izm, JAt,Izm0,Idx

  REAL(KREAL) :: sDihe, Pi, Sang

  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_mixed_frag")
  if (debug) Call Header("Entering Calc_mixed_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(FCart(Na),AtCart(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+1),DejaFait(Na),FCaf(na),FrozAt(na))

  Ind_Zmat=0

  ! Putting cart+frozen atoms into  cart array
  FCart=.FALSE.
  NCart=0
  Idx=1
  DO WHILE (Frozen(Idx).NE.0)
     If (Frozen(Idx).LT.0) THEN
        DO I=Frozen(Idx-1),abs(Frozen(Idx))
           IF (.NOT.FCart(I)) THEN
              NCart=NCart+1
              AtCart(NCart)=I
              FCart(I)=.TRUE.
           END IF
        END DO
     ELSEIF (.NOT.FCart(Frozen(Idx))) THEN
        FCart(Frozen(Idx))=.TRUE.
        NCart=NCart+1
        AtCart(NCart)=Frozen(Idx)
     END IF
  Idx=Idx+1
  END DO
  NFroz=NCart
  Idx=1
  DO WHILE (Cart(Idx).NE.0)
     If (Cart(Idx).LT.0) THEN
        DO I=Cart(Idx-1),abs(Cart(Idx))
           IF (.NOT.FCart(I)) THEN
              NCart=NCart+1
              AtCart(NCart)=I
              FCart(I)=.TRUE.
           END IF
        END DO
     ELSEIF (.NOT.FCart(Cart(Idx))) THEN
        FCart(Cart(Idx))=.TRUE.
        NCart=NCart+1
        AtCart(NCart)=Cart(Idx)
     END IF
  Idx=Idx+1
  END DO

  if (debug) THEN
     WRITE(*,'(1X,A,I4,A,I4,A)') "Found ",NFroz," frozen atoms, and a total of ",NCart," atoms described in cartesian coordinates"
     WRITE(*,*) "AtCart:",AtCart(1:NCart)
  END IF

  DejaFait=.FALSE.

  ! We cheat a lot: we will use ind_zmat et val_zmat to store
  ! the cartesian coordinates of the NCart atoms :-p
  DO I=1,NCart
     Iat=AtCart(I)
     ind_zmat(I,1)=Iat
     ind_zmat(I,2:5)=-1
     val_zmat(I,1)=X(Iat)
     val_zmat(I,2)=Y(Iat)
     val_zmat(I,3)=Z(Iat)
     DejaFait(Iat)=.TRUE.
     idx_zmat(Iat)=I
  END DO


  izm=NCart

  ! We now calculate the connections
  Liaisons=0
  LiaisonsBis=0

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

  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

  ! We get rid of bonds between cart atoms
  ! the trick is to do this on liaisons while keeping LiaisonsBis ok.

  LiaisonsIni=Liaisons
  LiaisonsBis=Liaisons

  DO I=1,NCart
     Iat=AtCart(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)
        IF (FCart(Jat)) THEN
           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 IF
     END DO
  END DO

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

  ! We 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
  Fragment=0
  NbAtFrag=0
  FrozBlock=0

  IdxFrag=0
  NbFrag=0

 Call Decomp_frag(na,liaisons,.NOT.Fcart,nbfrag,Fragment,NbAtFrag,FragAt)

 ! We compute FrozBlock now
  DO IFrag=1,NbFrag
     DO I=1,NbAtFrag(IFrag)
        Iat=FragAt(IFrag,I)
        IF (FCart(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," fragment(s) with cart atoms, out of",NbFrag," fragment(s)"


  if (debug) THEN
     WRITE(*,*) "Connections before adding fragments"
     DO IL=1,na
        IF (Liaisons(IL,0).NE.0)  WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
     END DO
  END IF

  ! We will now add the fragments containing non cart 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
  IZm=NCart
  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
        IaFaire=1
        FCaf(Iat)=.TRUE.
        DO WHILE (CaFaire(IaFaire).NE.0)
           n1=CaFaire(IaFaire)
           I1=idx_zmat(n1)
           n2=ind_zmat(I1,2)
           n3=ind_zmat(I1,3)

        if (debug) WRITE(*,1003) n1,(LIAISONS(n1,JL),JL=0,NMaxL)
           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
              if (debug) WRITE(*,1003) Iat,(LIAISONS(Iat,JL),JL=0,NMaxL)
              Call Annul(Liaisons,Iat,n1)
              if (debug) WRITE(*,1003) Iat,(LIAISONS(Iat,JL),JL=0,NMaxL)
              Liaisons(iat,0)=Liaisons(Iat,0)-1
           if (debug) WRITE(*,*) 'Iafaire,n1,DejaFait,idxzmat,n2,n3',IaFaire,n1, &
                DejaFait(n1),idx_zmat(n1),n2,3
           if (debug) WRITE(*,*) 'Cafaire - 0', &
                (CaFaire(J),J=Iafaire,IdxCAfaire)


              ! <- PFL 2007.Apr.16
              IF (.NOT.DejaFait(Iat)) THEN
                 ! we add it to the Zmat
!                 WRITE(*,*) "coucou"
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
!                 WRITE(*,*) "coucou 2"
                 izm=izm+1
                 IF (izm.EQ.2) THEN
!                 WRITE(*,*) "coucou 3"
                    ind_zmat(izm,1)=IAt
                    ind_zmat(izm,2)=n1
                    ind_zmat(izm,3:5)=0
                    if (n2.EQ.-1) n2=Iat
                 END IF

!                 WRITE(*,*) "coucou 4"
                 if ((izm.GE.3).AND.(n2.EQ.-1)) THEN
!                 WRITE(*,*) "coucou 5"
!                 WRITE(*,*) "coucou 3 bis"
                       ! This atom is defined in cartesian coordinates
                       ! and has not been used to put a non cart atom.
                       ! we will find the closest atom linked to this one amongst the atoms
                       ! already stored in ind_zmat
                       d=1e9
                       JAt=-1
                       DO I=1,NCart
                          If (ind_zmat(I,1).NE.n1) THEN
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
                             ! we take only atoms that are not aligned
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
                                Jat=ind_zmat(I,1)
                                d=norm2
                             END IF
                          END IF
                       END DO
                       if (JAt.EQ.-1) THEN
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
                          STOP
                       END IF
                       n2=JAt
                 END IF! izm>=3 and n2.eq.-1

                 IF (izm.eq.3) THEN
                    ind_zmat(izm,1)=Iat
                    ind_zmat(izm,2)=n1
                    ind_zmat(izm,3)=n2
                 END IF  ! izm=3

                 IF (izm.GE.4) THEN
!                                     WRITE(*,*) "coucou 5 bis"
                    IF (n3.EQ.-1) THEN
                       ! This atom is defined in cartesian coordinates
                       ! and has not been used to put a non cart atom.
                       ! we will find the closest atom linked to this one amongst the atoms
                       ! already stored in ind_zmat
                       call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
                       d=1e9
                       JAt=-1
                       DO I=1,NCart
                          If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
                             SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
                             if (debug) WRITE(*,*) "1-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
                             ! we take only atoms that are not aligned
                             if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
                                Jat=ind_zmat(I,1)
                                d=norm3
                             END IF
                          END IF
                       END DO
                       if (JAt.EQ.-1) THEN
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
                          STOP
                       END IF
                       n3=JAt
                    END IF
                    !                    ind_zmat(I1,3)=n3
                    Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
                 END IF ! izm>=4
!                 WRITE(*,*) "coucou 6"
           if (debug) THEN
              WRITE(*,*) 'ind_zmat 0',izm
              DO Ip=1,NCart
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(ip,1:4)
                 END DO
              SELECT CASE (NCart)
                 CASE (1)
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
                    if (izm.GE.3) &
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
                    Ibeg=4
                 CASE (2)
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
                    Ibeg=4
                 CASE DEFAULT
                    IBeg=NCart+1
                 END SELECT
              DO Ip=IBeg,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

                 idx_zmat(iat)=izm
!                 Call Annul(Liaisons,iat,n1)
!                 Liaisons(iat,0)=Liaisons(Iat,0)-1
                 !                  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)
           I1=idx_zmat(n1)
           n2=ind_zmat(I1,2)
           n3=ind_zmat(I1,3)
           if (debug) WRITE(*,*) 'Iafaire,n1,DejaFait,idxzmat,n2,n3',IaFaire,n1, &
                DejaFait(n1),idx_zmat(n1),n2,3
           if (debug) WRITE(*,*) 'Cafaire', &
                (CaFaire(J),J=Iafaire,IdxCAfaire)


           DO I3=1,Liaisons(n1,0)
              if (debug) WRITE(*,*) "Here, I3=",I3
              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
              ! <- PFL 2007.Apr.16
              IF (.NOT.DejaFait(Iat)) THEN
                 ! we add it to the Zmat
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
                 izm=izm+1
                 IF (izm.EQ.2) THEN
                    ind_zmat(izm,1)=IAt
                    ind_zmat(izm,2)=n1
                    ind_zmat(izm,3:5)=0
                 ELSE
                    IF (n2.EQ.-1) THEN
                       ! This atom is defined in cartesian coordinates
                       ! and has not been used to put a non cart atom.
                       ! we will find the closest atom linked to this one amongst the atoms
                       ! already stored in ind_zmat
                       d=1e9
                       JAt=-1
                       DO I=1,NCart
                          If (ind_zmat(I,1).NE.n1) THEN
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
                             ! we take only atoms that are not aligned
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
                                Jat=ind_zmat(I,1)
                                d=norm2
                             END IF
                          END IF
                       END DO
                       if (JAt.EQ.-1) THEN
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
                          STOP
                       END IF
                       n2=JAt
                    END IF
                    !                    ind_zmat(I1,2)=n2
                    if (izm.EQ.3) THEN
                       ind_zmat(izm,1)=Iat
                       ind_zmat(izm,2)=n1
                       ind_zmat(izm,3)=n2
                    ELSE
                       IF (n3.EQ.-1) THEN
                          ! This atom is defined in cartesian coordinates
                          ! and has not been used to put a non cart atom.
                          ! we will find the closest atom linked to this one amongst the atoms
                          ! already stored in ind_zmat
                          call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
                          d=1e9
                          JAt=-1
                          DO I=1,NCart
                             If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
                                call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
                                SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
                             if (debug) WRITE(*,*) "2-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3

                                ! we take only atoms that are not aligned
                                if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
                                   Jat=ind_zmat(I,1)
                                   d=norm3
                                END IF
                             END IF
                          END DO
                          if (JAt.EQ.-1) THEN
                             WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
                             STOP
                          END IF
                          n3=JAt
                       END IF
                       !                    ind_zmat(I1,3)=n3
                       Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
                    END IF
                 END IF
                 idx_zmat(iat)=izm
              !                  Call Annul(Liaisons,n1,iat)
                 n3=iat
                 DejaFait(Iat)=.TRUE.
                 if (debug) WRITE(*,*) "For no reason, I3=",I3
              END IF
           END DO
           IaFaire=IaFaire+1

           if (debug) THEN
              WRITE(*,*) 'ind_zmat 6',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(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
        I1=idx_zmat(n1)
        n2=ind_zmat(I1,2)
        n3=ind_zmat(I1,3)

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

        
        call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
        izm=izm+1
        IF (izm.EQ.2) THEN
           ind_zmat(izm,1)=IAt
           ind_zmat(izm,2)=n1
           ind_zmat(izm,3:5)=0
        ELSE
           IF (n2.EQ.-1) THEN
! This atom is defined in cartesian coordinates
! and has not been used to put a non cart atom.
! we will find the closest atom linked to this one amongst the atoms
! already stored in ind_zmat
              d=1e9
              JAt=-1
              DO I=1,NCart
                 If (ind_zmat(I,1).NE.n1) THEN
                    call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
                    SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
! we take only atoms that are not aligned
                    if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
                       Jat=ind_zmat(I,1)
                       d=norm2
                    END IF
                 END IF
              END DO
              if (JAt.EQ.-1) THEN
                 WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
                 STOP
              END IF
              n2=JAt
           END IF
!                    ind_zmat(I1,2)=n2
        END IF
        if (izm.EQ.3) THEN
           ind_zmat(izm,1)=Iat
           ind_zmat(izm,2)=n1
           ind_zmat(izm,3)=n2
        ELSE
           IF (n3.EQ.-1) THEN
! This atom is defined in cartesian coordinates
! and has not been used to put a non cart atom.
! we will find the closest atom linked to this one amongst the atoms
! already stored in ind_zmat
              call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
              d=1e9
              JAt=-1
              DO I=1,NCart
                 If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
                    call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
                    SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
                             if (debug) WRITE(*,*) "3-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3

! we take only atoms that are not aligned
                    if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
                       Jat=ind_zmat(I,1)
                       d=norm3
                    END IF
                 END IF
              END DO
              if (JAt.EQ.-1) THEN
                 WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
                 STOP
              END IF
              n3=JAt
           END IF
           !                    ind_zmat(I1,3)=n3
           Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
        END IF
        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)
           I1=idx_zmat(n1)
           n2=ind_zmat(I1,2)
           n3=ind_zmat(I1,3)

           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
              ! <- PFL 2007.Apr.16
              IF (.NOT.DejaFait(Iat)) THEN
                 ! we add it to the Zmat
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
                 izm=izm+1
                 IF (izm.EQ.2) THEN
                    ind_zmat(izm,1)=IAt
                    ind_zmat(izm,2)=n1
                    ind_zmat(izm,3:5)=0
                 ELSE
                    IF (n2.EQ.-1) THEN
! This atom is defined in cartesian coordinates
! and has not been used to put a non cart atom.
! we will find the closest atom linked to this one amongst the atoms
! already stored in ind_zmat
                       d=1e9
                       JAt=-1
                       DO I=1,NCart
                          If (ind_zmat(I,1).NE.n1) THEN
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
! we take only atoms that are not aligned
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
                                Jat=ind_zmat(I,1)
                                d=norm2
                             END IF
                          END IF
                       END DO
                       if (JAt.EQ.-1) THEN
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
                          STOP
                       END IF
                       n2=JAt
                    END IF
!                    ind_zmat(I1,2)=n2
                 END IF
                 if (izm.EQ.3) THEN
                    ind_zmat(izm,1)=Iat
                    ind_zmat(izm,2)=n1
                    ind_zmat(izm,3)=n2
                 ELSE
                    IF (n3.EQ.-1) THEN
! This atom is defined in cartesian coordinates
! and has not been used to put a non cart atom.
! we will find the closest atom linked to this one amongst the atoms
! already stored in ind_zmat
                       call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
                       d=1e9
                       JAt=-1
                       DO I=1,NCart
                          If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
                             SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
                             if (debug) WRITE(*,*) "4-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3

! we take only atoms that are not aligned
                             if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
                                Jat=ind_zmat(I,1)
                                d=norm3
                             END IF
                          END IF
                       END DO
                       if (JAt.EQ.-1) THEN
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
                          STOP
                       END IF
                       n3=JAt
                    END IF
!                    ind_zmat(I1,3)=n3
                    Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
                 END IF
                 idx_zmat(iat)=izm
                 n3=iat
                 DejaFait(Iat)=.TRUE.
              END IF
           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"

  SELECT CASE (NCart)
     CASE (1)
        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
     IBeg=4

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

     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
     IBeg=4
  CASE DEFAULT
     Ibeg=NCart+1
  END SELECT


  DO i=IBeg,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,NCart
        WRITE(*,'(1X,I5,3(1X,F11.6))') ind_zmat(i,1),(val_zmat(i,j),j=1,3)
     END DO
     DO I=NCart+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(*,*) "Deallocate FCart,AtCart,FCaf"
  DEALLOCATE(FCart,AtCart,FCaf)

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

END SUBROUTINE Calc_mixed_frag

