Chimie Théorique » Opt'n Path

      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

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

!  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 : Nat, NMaxL,  KINT, KREAL,Pi

      use Io_module

      IMPLICIT NONE

      INTEGER(KINT) :: Izm,I,IAf,IAt,IAt0,J,Jat

      INTEGER(KINT) :: n1, NMaxCF, NMinCP, VCF

      INTEGER(KINT) :: IAta,IAtd,IAtDh,IAti,IAtv,ICat

      INTEGER(KINT) :: ICP0, IdAt0, IdAtl

      INTEGER(KINT) :: Idx, IIat, IIatDh, IndAFaire, IndFin

      INTEGER(KINT) :: IdxAt0

      INTEGER(KINT) :: ITmp, IZm0, IZm1, IZm2, IZm3, IZm4, IzmTmp

      integer(KINT) :: na, atome(NAt), ind_zmat(NAt, 5)

      real(KREAL) ::  x(NAt), y(NAt), z(NAt)

      real(KREAL) ::  val_zmat(NAt,3)

!     ListAt contains the indices of frozen atoms

      LOGICAL ListAT(NAt), FFirst

      INTEGER(KINT) :: Natreat

      real(KREAL) ::  dist

      INTEGER(KINT) :: LIAISONS(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) :: NbCycle

      real(KREAL) ::  vx1,vy1,vz1,norm1

      real(KREAL) ::  vx2,vy2,vz2,norm2

      real(KREAL) ::  val

      Logical Debug, Done

      Logical DejaFait(NAt)

      INTERFACE

         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

      END INTERFACE

!     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