Chimie Théorique » Opt'n Path

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

  !

!----------------------------------------------------------------------

!  Copyright 2003-2014 Ecole Normale Supérieure de Lyon,

!  Centre National de la Recherche Scientifique,

!  Université Claude Bernard Lyon 1. All rights reserved.

!

!  This work is registered with the Agency for the Protection of Programs

!  as IDDN.FR.001.100009.000.S.P.2014.000.30625

!

!  Authors: P. Fleurat-Lessard, P. Dayal

!  Contact: optnpath@gmail.com

!

! This file is part of "Opt'n Path".

!

!  "Opt'n Path" is free software: you can redistribute it and/or modify

!  it under the terms of the GNU Affero General Public License as

!  published by the Free Software Foundation, either version 3 of the License,

!  or (at your option) any later version.

!

!  "Opt'n Path" is distributed in the hope that it will be useful,

!  but WITHOUT ANY WARRANTY; without even the implied warranty of

!

!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the

!  GNU Affero General Public License for more details.

!

!  You should have received a copy of the GNU Affero General Public License

!  along with "Opt'n Path". If not, see <http://www.gnu.org/licenses/>.

!

! Contact The Office of Technology Licensing, valorisation@ens-lyon.fr,

! for commercial licensing opportunities.

!----------------------------------------------------------------------

  Use Path_module, only : max_Z, NMaxL, Nom

  Use Io_module

  IMPLICIT NONE

! Parameters of the subroutine

! na: number of atoms in the system

  integer(KINT) ::  na

! atome: Mass number of the atoms of the system

  integer(KINT) ::  atome(na)

! ind_zmat: for "zmat" atoms contains the indices of reference atoms

  integer(KINT) ::  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)

  CHARACTER(5) :: AtName

  !     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 :: 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, d

  Logical Debug,DebugGaussian

  LOGICAL, ALLOCATABLE :: DejaFait(:),FCaf(:) !(na)

  Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen

  INTEGER(KINT) :: I, J, n1, n2, n3, n4, IAt, IL, JL, IFrag, ITmp

  INTEGER(KINT) :: I3, I1, Ip, IFragFroz, IBeg

  INTEGER(KINT) :: I0, Izm, JAt,Izm0,Idx

  REAL(KREAL) :: 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")

  debugGaussian=valid("zmat_gaussian")

  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,vx2,vy2,vz2, &

          vx4,vy4,vz4,norm4)

     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &

          vx5,vy5,vz5,norm5)

     val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &

          vx2,vy2,vz2,norm2)

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

!11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))