Statistiques
| Révision :

root / src / ZmatBuild.f90 @ 1

Historique | Voir | Annoter | Télécharger (31,29 ko)

1 1 pfleura2
      SUBROUTINE ZmatBuild(na,atome,ind_zmat,val_zmat,x,y,z,izm, &
2 1 pfleura2
           Liaisons, LIes_CP,LIes_CF,ListAt,DejaFait,debug)
3 1 pfleura2
4 1 pfleura2
!
5 1 pfleura2
!     This second version enables the user to freeze some atoms
6 1 pfleura2
!     the frozen atoms indices are listed in the frozen array.
7 1 pfleura2
!
8 1 pfleura2
!     The idea is surely too stupid to work all the time... but here it is
9 1 pfleura2
!     we first construct the zmat using only the frozen atoms, and then
10 1 pfleura2
!     we add the other atoms on top of the first ones...
11 1 pfleura2
!
12 1 pfleura2
! 12.06.06
13 1 pfleura2
! Small modification: there was an inconsistency: izm was increased
14 1 pfleura2
! after the atom was added to ind_zmat for the three first atoms,
15 1 pfleura2
! but before for all others. Now, it starts at zero and it
16 1 pfleura2
! is increased before the atom is added to ind_zmat
17 1 pfleura2
!     IMPLICIT NONE
18 1 pfleura2
19 1 pfleura2
      use Path_module, only : Nat, NMaxL, max_Z, KINT, KREAL
20 1 pfleura2
21 1 pfleura2
22 1 pfleura2
      INTEGER(KINT) :: Izm
23 1 pfleura2
      integer(KINT) :: na,atome(NAt),at,ind_zmat(NAt,5),long
24 1 pfleura2
      real(KREAL) ::  x(NAt),y(NAt),z(NAt),fact
25 1 pfleura2
      real(KREAL) ::  val_zmat(NAt,3)
26 1 pfleura2
!     ListAt contains the indices of frozen atoms
27 1 pfleura2
      LOGICAL ListAT(NAt),FIzm1,FFirst
28 1 pfleura2
      INTEGER(KINT) :: Natreat
29 1 pfleura2
30 1 pfleura2
      real(KREAL) ::  vx,vy,vz,dist
31 1 pfleura2
      INTEGER(KINT) :: LIAISONS(NAt,0:NMaxL),Nbli,Nblj
32 1 pfleura2
      INTEGER(KINT) :: LiaisonsBis(NAt,0:NMaxL)
33 1 pfleura2
      INTEGER(KINT) :: CAFaire(NAt)
34 1 pfleura2
      INTEGER(KINT) :: LieS_CP(NAt,0:NMaxL),LieS_CF(NAt,0:NMaxL)
35 1 pfleura2
      INTEGER(KINT) :: AtCP0(NAt),NbAt0,NbAt0r
36 1 pfleura2
      Integer(KINT) :: AtTypCycl(Max_Z), NbAt3(NAt,2),CyCleAt(NAt)
37 1 pfleura2
      Integer(KINT) :: NbCycle
38 1 pfleura2
      real(KREAL) ::  vx1,vy1,vz1,norm1
39 1 pfleura2
      real(KREAL) ::  vx2,vy2,vz2,norm2
40 1 pfleura2
      real(KREAL) ::  vx3,vy3,vz3,norm3
41 1 pfleura2
      real(KREAL) ::  vx4,vy4,vz4,norm4
42 1 pfleura2
      real(KREAL) ::  vx5,vy5,vz5,norm5
43 1 pfleura2
      real(KREAL) ::  val,val_d
44 1 pfleura2
      Logical AtTerm,Debug,Done
45 1 pfleura2
      Logical DejaFait(NAt)
46 1 pfleura2
      LOGICAL F1213, F1223,F1323,FTmp
47 1 pfleura2
48 1 pfleura2
!     As we may have to treat only partial molecules, it may happen
49 1 pfleura2
!     that there are no central atoms... so we first check this by looking
50 1 pfleura2
!     for the least number of CP atoms
51 1 pfleura2
      FFirst=.TRUE.
52 1 pfleura2
      NMinCP=na
53 1 pfleura2
      NMaxCF=-1
54 1 pfleura2
      NAtreat=0
55 1 pfleura2
      Izm0=Izm+1
56 1 pfleura2
57 1 pfleura2
      DO I=1,na
58 1 pfleura2
         if (DEBUG) WRITE(*,*) "DBG ZmatBuild i,listat",i,listAt(i)
59 1 pfleura2
         IF (ListAt(I).AND.(.NOT.(DejaFait(I)))) THEN
60 1 pfleura2
            IF (Lies_CP(I,0).lt.NMinCP)      NMinCP=Lies_CP(I,0)
61 1 pfleura2
            IF (Lies_CF(I,0).gt.NMaxCF)      NMaxCF=Lies_CF(I,0)
62 1 pfleura2
            NATreat=NATreat+1
63 1 pfleura2
         END IF
64 1 pfleura2
      END DO
65 1 pfleura2
      if (debug) WRITE(*,*) 'NMinCP=',NMinCP
66 1 pfleura2
      if (debug) WRITE(*,*) 'NMaxCF=',NMaxCF
67 1 pfleura2
      if (debug) WRITE(*,*) 'NATreat=',NAtreat
68 1 pfleura2
69 1 pfleura2
70 1 pfleura2
      IF (Debug) THEN
71 1 pfleura2
         WRITE(*,*) '------------------ Zmat Build -------------------'
72 1 pfleura2
         WRITE(*,*) 'LIAISONS'
73 1 pfleura2
         DO I=1,na
74 1 pfleura2
            WRITE(*,1003) i,(LIAISONS(I,J),J=0,NMaxL)
75 1 pfleura2
         END DO
76 1 pfleura2
77 1 pfleura2
         WRITE(*,*) 'LIes_CF'
78 1 pfleura2
         DO I=1,na
79 1 pfleura2
            WRITE(*,1003) i,(LIeS_CF(I,J),J=0,NMaxL)
80 1 pfleura2
         END DO
81 1 pfleura2
82 1 pfleura2
         WRITE(*,*) 'LIes_CP'
83 1 pfleura2
         DO I=1,na
84 1 pfleura2
            WRITE(*,1003) i,(LIeS_CP(I,J),J=0,NMaxL)
85 1 pfleura2
         END DO
86 1 pfleura2
         WRITE(*,*) '-- Zmat Build : only _To treat_ atoms------------'
87 1 pfleura2
         WRITE(*,*) 'LIAISONS'
88 1 pfleura2
         DO I=1,na
89 1 pfleura2
            IF (ListAt(I).AND.(.NOT.DejaFait(I))) &
90 1 pfleura2
                 WRITE(*,1003) i,(LIAISONS(I,J),J=0,NMaxL)
91 1 pfleura2
         END DO
92 1 pfleura2
93 1 pfleura2
         WRITE(*,*) 'LIes_CF'
94 1 pfleura2
         DO I=1,na
95 1 pfleura2
            IF (ListAt(I).AND.(.NOT.DejaFait(I))) &
96 1 pfleura2
                WRITE(*,1003) i,(LIeS_CF(I,J),J=0,NMaxL)
97 1 pfleura2
         END DO
98 1 pfleura2
99 1 pfleura2
         WRITE(*,*) 'LIes_CP'
100 1 pfleura2
         DO I=1,na
101 1 pfleura2
            IF (ListAt(I).AND.(.NOT.DejaFait(I))) &
102 1 pfleura2
                WRITE(*,1003) i,(LIeS_CP(I,J),J=0,NMaxL)
103 1 pfleura2
         END DO
104 1 pfleura2
105 1 pfleura2
      END IF
106 1 pfleura2
107 1 pfleura2
108 1 pfleura2
      if (NAtreat.EQ.0) THEN
109 1 pfleura2
         WRITE(*,*) "Foutage de gueule : NAtTreat=0"
110 1 pfleura2
         RETURN
111 1 pfleura2
      END IF
112 1 pfleura2
113 1 pfleura2
      if (debug) THEN
114 1 pfleura2
         WRITE(*,*) 'DBG ZmatBuil, izm=',izm
115 1 pfleura2
         DO I=1,na
116 1 pfleura2
            WRITE(*,'(1X,I5,1X,L4,1X,L4,12(1X,I4))') i,ListAt(I) &
117 1 pfleura2
                 ,DejaFait(I), &
118 1 pfleura2
                 (Liaisons(I,J),J=0,Liaisons(I,0))
119 1 pfleura2
         END DO
120 1 pfleura2
      END IF
121 1 pfleura2
122 1 pfleura2
!!!   atraiter NMaxCF=0 :que des atomes separes...
123 1 pfleura2
!     mais faut-il reellement faire un cas a part ?
124 1 pfleura2
125 1 pfleura2
!     On compte le nb d'atomes sans atomes CP (centripetes)
126 1 pfleura2
      NbAt0=0
127 1 pfleura2
      DO I=1,na
128 1 pfleura2
         IF (ListAt(I).AND.(.NOT.DejaFait(I)).AND. &
129 1 pfleura2
              (Lies_CP(I,0).eq.NMinCP)) THEN
130 1 pfleura2
            NbAt0=NbAt0+1
131 1 pfleura2
            AtCP0(NbAt0)=I
132 1 pfleura2
            IF (DEBUG) WRITE(*,*) "Center atom",NbAt0,I,Atome(I)
133 1 pfleura2
         END IF
134 1 pfleura2
      END DO
135 1 pfleura2
      AtCP0(NbAt0+1)=0
136 1 pfleura2
      NbAt0r=NbAt0
137 1 pfleura2
138 1 pfleura2
      IF (Debug)  WRITE(*,*) 'DBG ZmatBuld - NCP0,AtCP0 ',NbAt0, &
139 1 pfleura2
           (AtCP0(i),i=1,NbAt0)
140 1 pfleura2
141 1 pfleura2
!     On va les traiter tous en construisant les molecules en partant d'eux
142 1 pfleura2
!     Le premier atome sans CP est different des autres car il fixe
143 1 pfleura2
!     l'origine
144 1 pfleura2
145 1 pfleura2
!First atom
146 1 pfleura2
      IF (Izm==0) THEN
147 1 pfleura2
         FFirst=.FALSE.
148 1 pfleura2
         if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=1'
149 1 pfleura2
!     IZm=1
150 1 pfleura2
!     Boucle pour trouver celui qui a le plus d'atomes CF
151 1 pfleura2
         IdAt0=0
152 1 pfleura2
         VCF=-1
153 1 pfleura2
         DO I=1,NbAt0
154 1 pfleura2
!     WRITE(*,*) 'DBG ZmatB',I,AtCP0(I)
155 1 pfleura2
            if (AtCP0(I).NE.0) THEN
156 1 pfleura2
               ICat=AtCP0(I)
157 1 pfleura2
               if (Lies_CF(ICAt,0).gt.VCF) THEN
158 1 pfleura2
                  IdAt0=ICat
159 1 pfleura2
                  IdxAt0=I
160 1 pfleura2
                  VCF=Lies_CF(ICAt,0)
161 1 pfleura2
               END IF
162 1 pfleura2
            END IF
163 1 pfleura2
         END DO
164 1 pfleura2
165 1 pfleura2
!     IF (debug) WRITE(*,*) 'DBG ZmatBuil - VCF, IdxAt0,IdAt0',
166 1 pfleura2
!     &        VCF, IdxAt0,IdAt0
167 1 pfleura2
!     all atoms with NMinCP CP links do not have CF links... not a good choice to start
168 1 pfleura2
!     building the molecule(s) around them... we increase NMinCP
169 1 pfleura2
         AtCP0(IdxAt0)=0
170 1 pfleura2
         NbAt0r=NbAt0r-1
171 1 pfleura2
172 1 pfleura2
!     WRITE(IOOUT,*) 'Atome 1', IdAt0, Nom(Atome(IdAt0))
173 1 pfleura2
         Izm1=IdAt0
174 1 pfleura2
         Izm=Izm+1
175 1 pfleura2
176 1 pfleura2
         ind_zmat(izm,1)=Izm1
177 1 pfleura2
         ind_zmat(izm,2)=0
178 1 pfleura2
         ind_zmat(izm,3)=0
179 1 pfleura2
         ind_zmat(izm,4)=0
180 1 pfleura2
         ind_zmat(izm,5)=0
181 1 pfleura2
         val_zmat(izm,1)=0.0
182 1 pfleura2
         val_zmat(izm,2)=0.0
183 1 pfleura2
         val_zmat(izm,3)=0.0
184 1 pfleura2
         DejaFait(Izm1)=.True.
185 1 pfleura2
186 1 pfleura2
      END IF                    ! end of izm==1 test
187 1 pfleura2
188 1 pfleura2
!     On construit la premiere couche qui est speciale car elle fixe les
189 1 pfleura2
!     axes.
190 1 pfleura2
!     Plusieurs cas sont possibles suivant le nb d'atomes CF
191 1 pfleura2
      IdAt0=ind_zmat(1,1)
192 1 pfleura2
      Izm1=IdAt0
193 1 pfleura2
194 1 pfleura2
!Second atom
195 1 pfleura2
      IF ((izm==1).AND.(NAtreat.ge.2))  THEN
196 1 pfleura2
         Done=.FALSE.
197 1 pfleura2
         if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=2'
198 1 pfleura2
!     Two cases: 1) We are constructing the first fragment: FIzm1=.T.
199 1 pfleura2
!     2) We are constructing the second fragment: FIzm1=.F.
200 1 pfleura2
!     in this case, we select one CP0 atom to start with
201 1 pfleura2
         IF (FFirst) THEN
202 1 pfleura2
            FFirst=.FALSE.
203 1 pfleura2
!     This is the first atom to be dealt with
204 1 pfleura2
!     We look for a CP0
205 1 pfleura2
            If (NbAt0r.ge.1) THEN
206 1 pfleura2
               IdAt0=0
207 1 pfleura2
               VCF=-1
208 1 pfleura2
               WRITE(*,*) 'DBG  ZBuild Izm=2:AtCP0',(AtCP0(I),I=1,NbAt0)
209 1 pfleura2
               DO I=1,NbAt0
210 1 pfleura2
                  if (AtCP0(I).NE.0) THEN
211 1 pfleura2
                     ICat=AtCP0(I)
212 1 pfleura2
                     if (Lies_CF(I,0).gt.VCF) THEN
213 1 pfleura2
                        Izm2=ICat
214 1 pfleura2
                        IdxAt0=I
215 1 pfleura2
                        VCF=Lies_CF(ICAt,0)
216 1 pfleura2
                     END IF
217 1 pfleura2
                  END IF
218 1 pfleura2
               END DO
219 1 pfleura2
               WRITE(*,*) 'DBG ZBuild Izm=2:izm2,IdxAt0,VCF,AtCP0(1)', &
220 1 pfleura2
                    izm2,IdxAt0,VCF,AtCP0(1)
221 1 pfleura2
               NbAt0r=NbAt0r-1
222 1 pfleura2
               Done=.TRUE.
223 1 pfleura2
            END IF
224 1 pfleura2
225 1 pfleura2
!     If we do not have a CP0 available the other tests are identical
226 1 pfleura2
!     for cases 1 and 2...
227 1 pfleura2
         END IF
228 1 pfleura2
         IF (.NOT.DONE) THEN
229 1 pfleura2
            IF (Lies_CF(IdAt0,0).ge.1) THEN
230 1 pfleura2
               IZm2=Lies_CF(IdAt0,1)
231 1 pfleura2
               WRITE(*,*) 'DBG ZBuild Lies_CF(IdAt0,0).ge.1', &
232 1 pfleura2
                    Lies_CF(IdAt0,0),izm2
233 1 pfleura2
            ELSE
234 1 pfleura2
!     First atom is not linked to anything... we look for the second CP0 atom
235 1 pfleura2
!     if we have one..
236 1 pfleura2
               If (NbAt0r.ge.1) THEN
237 1 pfleura2
                  IdAt0=0
238 1 pfleura2
                  VCF=-1
239 1 pfleura2
                  IF (DEBUG) WRITE(*,*) 'DBG  ZBuild Izm=2:AtCP0', &
240 1 pfleura2
                       (AtCP0(I),I=1,NbAt0)
241 1 pfleura2
                  DO I=1,NbAt0
242 1 pfleura2
                     if (AtCP0(I).NE.0) THEN
243 1 pfleura2
                        ICat=AtCP0(I)
244 1 pfleura2
                        if (Lies_CF(I,0).gt.VCF) THEN
245 1 pfleura2
                           Izm2=ICat
246 1 pfleura2
                           IdxAt0=I
247 1 pfleura2
                           VCF=Lies_CF(ICAt,0)
248 1 pfleura2
                        END IF
249 1 pfleura2
                     END IF
250 1 pfleura2
                  END DO
251 1 pfleura2
                  IF (DEBUG) WRITE(*,*) 'DBG ZBuild Izm=2:izm2,IdxAt0,VCF,AtCP0(1)', &
252 1 pfleura2
                       izm2,IdxAt0,VCF,AtCP0(1)
253 1 pfleura2
                  NbAt0r=NbAt0r-1
254 1 pfleura2
               ELSE
255 1 pfleura2
!     we do not have another CP0 atom... but we know that there
256 1 pfleura2
!     is at least two atoms... we look for the closest atom to Izm1
257 1 pfleura2
                  Izm2=0
258 1 pfleura2
                  Dist=1.D99
259 1 pfleura2
                  DO I=1,Na
260 1 pfleura2
                     IF (ListAt(I).AND.(.NOT.Dejafait(I))) THEN
261 1 pfleura2
                        CALL vecteur(IZm1,I,x,y,z,vx1,vy1,vz1,norm1)
262 1 pfleura2
                        if (Norm1.lt.Dist) IZm2=I
263 1 pfleura2
                     END IF
264 1 pfleura2
                  END DO
265 1 pfleura2
               END IF
266 1 pfleura2
            END IF
267 1 pfleura2
         END IF
268 1 pfleura2
         CALL vecteur(IZm1,IZm2,x,y,z,vx1,vy1,vz1,norm1)
269 1 pfleura2
270 1 pfleura2
         Izm=Izm+1
271 1 pfleura2
272 1 pfleura2
         ind_zmat(izm,1)=izm2
273 1 pfleura2
         ind_zmat(izm,2)=izm1
274 1 pfleura2
         ind_zmat(izm,3)=0
275 1 pfleura2
         ind_zmat(izm,4)=0
276 1 pfleura2
         ind_zmat(izm,5)=0
277 1 pfleura2
         val_zmat(izm,1)=norm1
278 1 pfleura2
         val_zmat(izm,2)=0.0
279 1 pfleura2
         val_zmat(izm,3)=0.0
280 1 pfleura2
         DejaFait(Izm2)=.TRUE.
281 1 pfleura2
282 1 pfleura2
      END IF                    ! end of izm==2 test
283 1 pfleura2
284 1 pfleura2
      Izm1=ind_zmat(1,1)
285 1 pfleura2
      Izm2=ind_zmat(2,1)
286 1 pfleura2
287 1 pfleura2
! Third atom
288 1 pfleura2
289 1 pfleura2
      IF ((Izm==2).AND.(NAtreat.ge.3)) THEN
290 1 pfleura2
         Done=.FALSE.
291 1 pfleura2
         if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=3',FFirst
292 1 pfleura2
!     Two cases: 1) We are constructing the first fragment: FIzm1=.T.
293 1 pfleura2
!     2) We are constructing the second fragment: FIzm1=.F.
294 1 pfleura2
!     in this case, we select one CP0 atom to start with
295 1 pfleura2
         IF (FFirst) THEN
296 1 pfleura2
            FFirst=.FALSE.
297 1 pfleura2
!     This is the first atom to be dealt with
298 1 pfleura2
!     We look for a CP0
299 1 pfleura2
            If (NbAt0r.ge.1) THEN
300 1 pfleura2
               IdAt0=0
301 1 pfleura2
               VCF=-1
302 1 pfleura2
               DO I=1,NbAt0
303 1 pfleura2
                  if (AtCP0(I).NE.0) THEN
304 1 pfleura2
                     ICat=AtCP0(I)
305 1 pfleura2
                     if (Lies_CF(I,0).gt.VCF) THEN
306 1 pfleura2
                        Izm3=ICat
307 1 pfleura2
                        IdxAt0=I
308 1 pfleura2
                        VCF=Lies_CF(ICAt,0)
309 1 pfleura2
                     END IF
310 1 pfleura2
                  END IF
311 1 pfleura2
               END DO
312 1 pfleura2
               NbAt0r=NbAt0r-1
313 1 pfleura2
               Done=.TRUE.
314 1 pfleura2
            END IF
315 1 pfleura2
!     If we do not have a CP0 available the other tests are identical
316 1 pfleura2
!     for cases 1 and 2...
317 1 pfleura2
         END IF
318 1 pfleura2
         IF (.NOT.Done) THEN
319 1 pfleura2
            WRITE(*,*) 'DBG ZmatBuild: Done false for izm=3'
320 1 pfleura2
            IF (Lies_CF(izm1,0).ge.2) THEN
321 1 pfleura2
!     easiest case: the first atom is linked to at least 2 atoms.
322 1 pfleura2
               WRITE(*,*) 'DBG ZmtB-izm3:Lies_CF(izm1,0)>=2 '
323 1 pfleura2
               I=1
324 1 pfleura2
               DO WHILE ((.NOT.ListAt(Lies_CF(izm1,I))) &
325 1 pfleura2
                    .OR.(DejaFait(Lies_CF(izm1,I))))
326 1 pfleura2
                  I=I+1
327 1 pfleura2
               END DO
328 1 pfleura2
               izm3=Lies_CF(Izm1,I)
329 1 pfleura2
               Done=.TRUE.
330 1 pfleura2
            ELSEIF (Lies_CF(Izm2,0).ge.1) THEN
331 1 pfleura2
!     a bit more complex: first atom is linked to one atom (Izm2)
332 1 pfleura2
!     but luckily 2nd atom is linked to at least one atom
333 1 pfleura2
               I=1
334 1 pfleura2
               DO WHILE ((.NOT.ListAt(Lies_CF(izm2,I))) &
335 1 pfleura2
                    .OR.(DejaFait(Lies_CF(izm2,I))).AND.(I.le.Lies_CF(Izm2,0)))
336 1 pfleura2
                  I=I+1
337 1 pfleura2
               END DO
338 1 pfleura2
               IF (I.LE.Lies_CF(Izm2,0)) THEN
339 1 pfleura2
                  Izm3=Lies_CF(izm2,I)
340 1 pfleura2
               WRITE(*,*) 'DBG ZmtB-izm3:Lies_CF(izm2,0)>=1 and ok  '
341 1 pfleura2
!     We exchange Izm1 and Izm2 because the logical order is Izm3 linked
342 1 pfleura2
!     to Izm2  and not to Izm1 as is the default later
343 1 pfleura2
                  IzmTmp=Izm2
344 1 pfleura2
                  Izm2=Izm1
345 1 pfleura2
                  Izm1=IzmTmp
346 1 pfleura2
                  Done=.TRUE.
347 1 pfleura2
               END IF
348 1 pfleura2
            END IF
349 1 pfleura2
         IF (.NOT.Done) THEN
350 1 pfleura2
!     None of the first two atoms has links to a third atom...
351 1 pfleura2
!     we take it from the CP0 if possible...
352 1 pfleura2
               If (NbAt0r.ge.1) THEN
353 1 pfleura2
                  WRITE(*,*) 'DBG ZmatB izm=3',NbAt0r
354 1 pfleura2
                  IdAt0=0
355 1 pfleura2
                  VCF=-1
356 1 pfleura2
                  DO I=1,NbAt0
357 1 pfleura2
                     if (AtCP0(I).NE.0) THEN
358 1 pfleura2
                        ICat=AtCP0(I)
359 1 pfleura2
                        if (Lies_CF(I,0).gt.VCF) THEN
360 1 pfleura2
                           Izm3=ICat
361 1 pfleura2
                           IdxAt0=I
362 1 pfleura2
                           VCF=Lies_CF(ICAt,0)
363 1 pfleura2
                        END IF
364 1 pfleura2
                     END IF
365 1 pfleura2
                  END DO
366 1 pfleura2
                  NbAt0r=NbAt0r-1
367 1 pfleura2
                  WRITE(*,*) 'DBG ZmatB izm=3',NbAt0r,Izm3
368 1 pfleura2
               ELSE
369 1 pfleura2
! Or the  atom closest to 1
370 1 pfleura2
                  Izm3=0
371 1 pfleura2
                  Dist=1.D99
372 1 pfleura2
                  DO I=1,Na
373 1 pfleura2
                     IF (ListAt(I).AND.(.NOT.Dejafait(I))) THEN
374 1 pfleura2
                        CALL vecteur(IZm1,I,x,y,z,vx1,vy1,vz1,norm1)
375 1 pfleura2
                        if (Norm1.lt.Dist) IZm3=I
376 1 pfleura2
                     END IF
377 1 pfleura2
                  END DO
378 1 pfleura2
               END IF
379 1 pfleura2
            END IF
380 1 pfleura2
         END IF
381 1 pfleura2
382 1 pfleura2
         CALL vecteur(IZm1,IZm2,x,y,z,vx1,vy1,vz1,norm1)
383 1 pfleura2
         CALL vecteur(izm1,izm3,x,y,z,vx2,vy2,vz2,norm2)
384 1 pfleura2
         val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
385 1 pfleura2
386 1 pfleura2
!     write(*,11) n1,n2,norm1,n3,val
387 1 pfleura2
388 1 pfleura2
         Izm=Izm+1
389 1 pfleura2
390 1 pfleura2
         ind_zmat(izm,1)=izm3
391 1 pfleura2
         ind_zmat(izm,2)=izm1
392 1 pfleura2
         ind_zmat(izm,3)=izm2
393 1 pfleura2
         ind_zmat(izm,4)=0
394 1 pfleura2
         ind_zmat(izm,5)=0
395 1 pfleura2
         val_zmat(izm,1)=norm2
396 1 pfleura2
         val_zmat(izm,2)=val
397 1 pfleura2
         val_zmat(izm,3)=0.0
398 1 pfleura2
         DejaFait(Izm3)=.TRUE.
399 1 pfleura2
400 1 pfleura2
      END IF                    ! end of test izm=3
401 1 pfleura2
402 1 pfleura2
!     We add all atoms CF atoms linked to the already present atoms in the zmat to
403 1 pfleura2
!     the "to do" list:
404 1 pfleura2
!     Les atomes CF lies a IdAt0 sont mis en attente pour
405 1 pfleura2
!     etre traites
406 1 pfleura2
407 1 pfleura2
!     For 'historical' reasons, the first atom links have to be dealt with
408 1 pfleura2
!     in a special way... now !
409 1 pfleura2
410 1 pfleura2
      if (Izm.ge.NaTreat) Return
411 1 pfleura2
412 1 pfleura2
413 1 pfleura2
      WRITE(*,*) 'DBG ZmatBuild Izm0',Izm0
414 1 pfleura2
415 1 pfleura2
      if (debug) THEN
416 1 pfleura2
         WRITE(*,*) 'ind_zmat'
417 1 pfleura2
         DO I=1,na
418 1 pfleura2
            WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
419 1 pfleura2
         END DO
420 1 pfleura2
      END IF
421 1 pfleura2
422 1 pfleura2
      IF (FFIrst) GOTO 9753
423 1 pfleura2
424 1 pfleura2
         I=Ind_zmat(Izm0,1)
425 1 pfleura2
         DO IdAtl=1,Lies_CF(I,0)
426 1 pfleura2
            Izm4= Lies_CF(I,IdAtl)
427 1 pfleura2
!      WRITE(IOOUT,*) 'Atome ',Izm+1,Izm4,Nom(Atome(Izm4)),Izm1,
428 1 pfleura2
!     $Izm2,Izm3
429 1 pfleura2
            IF (ListAt(Izm4).AND.(.NOT.DejaFait(Izm4))) THEN
430 1 pfleura2
               Izm=Izm+1
431 1 pfleura2
               Call AddLigne(Izm,Izm4,Izm1,Izm2,Izm3,ind_zmat,val_zmat, &
432 1 pfleura2
                 x,y,z)
433 1 pfleura2
               DejaFait(Izm4)=.TRUE.
434 1 pfleura2
            END IF
435 1 pfleura2
         END DO
436 1 pfleura2
!      END IF
437 1 pfleura2
438 1 pfleura2
      if (debug) THEN
439 1 pfleura2
         WRITE(*,*) 'ind_zmat'
440 1 pfleura2
         DO I=1,na
441 1 pfleura2
            WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
442 1 pfleura2
         END DO
443 1 pfleura2
       END IF
444 1 pfleura2
445 1 pfleura2
      if (Izm.ge.NaTreat) Return
446 1 pfleura2
447 1 pfleura2
      IndFin=1
448 1 pfleura2
      DO I=Izm0,Izm
449 1 pfleura2
         Iat=ind_zmat(I,1)
450 1 pfleura2
!         Do iaf=1,Lies_CF(Iat,0)
451 1 pfleura2
!!     if (ListAt(Lies_CF(IAt,iaf)).AND.
452 1 pfleura2
!!     &              (.NOT.DejaFait(Lies_CF(IAt,iaf))))  THEN
453 1 pfleura2
!            if (ListAt(Lies_CF(IAt,iaf))) THEN
454 1 pfleura2
!               CAfaire(IndFin)=Lies_CF(IAt,Iaf)
455 1 pfleura2
               CaFaire(IndFin)=Iat
456 1 pfleura2
               IndFin=IndFin+1
457 1 pfleura2
!            END IF
458 1 pfleura2
!         END DO
459 1 pfleura2
      END DO
460 1 pfleura2
      CAfaire(IndFin)=0
461 1 pfleura2
462 1 pfleura2
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
463 1 pfleura2
464 1 pfleura2
      if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=4'
465 1 pfleura2
466 1 pfleura2
!     on a cree la premiere couche autour du premier centre
467 1 pfleura2
!     reste a finir les autres couches
468 1 pfleura2
      IndAFaire=1
469 1 pfleura2
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
470 1 pfleura2
      Do WHILE (Cafaire(IndAFaire).ne.0)
471 1 pfleura2
         IAt0=Cafaire(IndAfaire)
472 1 pfleura2
         if (debug) WRITE(*,*) "IndAFaire,IAt0,Lies_CF(IAt0,0)", &
473 1 pfleura2
              IndAFaire,IAt0,Lies_CF(IAt0,0)
474 1 pfleura2
         if (Lies_CF(IAt0,0).ge.1) THEN
475 1 pfleura2
!     IAt0 est l'atome pour lequel on construit la couche suivante
476 1 pfleura2
!     le premier atome lie est particulier car il definit l'orientation
477 1 pfleura2
!     de ce fragment
478 1 pfleura2
            IAti=Lies_CF(IAt0,1)
479 1 pfleura2
      WRITE(IOOUT,*) 'IAti:',IAti
480 1 pfleura2
            IAtd=IAt0
481 1 pfleura2
      WRITE(IOOUT,*) 'IAtd:',IAtd
482 1 pfleura2
            IAtv=Lies_CP(IAt0,1)
483 1 pfleura2
      WRITE(IOOUT,*) 'IAtv:',IAtv
484 1 pfleura2
            IIAtdh=1
485 1 pfleura2
            IAtdh=Lies_CF(IAtv,IIatdh)
486 1 pfleura2
            DO WHILE (Iatdh.eq.IAt0)
487 1 pfleura2
               IIatdh=IIatdh+1
488 1 pfleura2
               IAtdh=Lies_CF(IAtv,IIatdh)
489 1 pfleura2
            END DO
490 1 pfleura2
            IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
491 1 pfleura2
      WRITE(IOOUT,*) 'IAtdh:',IAtdh
492 1 pfleura2
            IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
493 1 pfleura2
               Izm=Izm+1
494 1 pfleura2
               if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
495 1 pfleura2
                    izm,IAti,IAtd,IAtv,IAtdh
496 1 pfleura2
               Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
497 1 pfleura2
                    ind_zmat,val_zmat ,x,y,z)
498 1 pfleura2
               DejaFait(IAti)=.TRUE.
499 1 pfleura2
            END IF
500 1 pfleura2
501 1 pfleura2
            CAfaire(IndFin)=Lies_CF(IAt0,1)
502 1 pfleura2
            IndFin=IndFin+1
503 1 pfleura2
!     Le traitement des autres est plus facile
504 1 pfleura2
            IAtdh=Lies_CF(IAt0,1)
505 1 pfleura2
            DO IAta=2,Lies_CF(IAt0,0)
506 1 pfleura2
               IAtI=Lies_CF(IAt0,IAta)
507 1 pfleura2
               CAfaire(IndFin)=IAtI
508 1 pfleura2
               IndFin=IndFin+1
509 1 pfleura2
510 1 pfleura2
               IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
511 1 pfleura2
                  Izm=Izm+1
512 1 pfleura2
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh,ind_zmat &
513 1 pfleura2
                       ,val_zmat,x,y,z)
514 1 pfleura2
                  DejaFait(IAti)=.TRUE.
515 1 pfleura2
               END IF
516 1 pfleura2
            END DO
517 1 pfleura2
            CAfaire(IndFin)=0
518 1 pfleura2
         END IF
519 1 pfleura2
         IndAFaire=IndAFaire+1
520 1 pfleura2
         if (debug) WRITE(IOOUT,*) "IndAfaire,CaFaire:",IndAfaire, &
521 1 pfleura2
            CaFaire
522 1 pfleura2
      END DO
523 1 pfleura2
524 1 pfleura2
      IndFin=1
525 1 pfleura2
      DO I=1,Izm0
526 1 pfleura2
         Iat=ind_zmat(I,1)
527 1 pfleura2
!         Do iaf=1,Lies_CF(Iat,0)
528 1 pfleura2
         if (ListAt(IAt).AND. &
529 1 pfleura2
                    (.NOT.DejaFait(IAt)))  THEN
530 1 pfleura2
!            if (ListAt(Lies_CF(IAt,iaf))) THEN
531 1 pfleura2
!               CAfaire(IndFin)=Lies_CF(IAt,Iaf)
532 1 pfleura2
               CaFaire(IndFin)=Iat
533 1 pfleura2
               IndFin=IndFin+1
534 1 pfleura2
            END IF
535 1 pfleura2
!         END DO
536 1 pfleura2
      END DO
537 1 pfleura2
      CAfaire(IndFin)=0
538 1 pfleura2
539 1 pfleura2
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
540 1 pfleura2
541 1 pfleura2
      if (debug) WRITE(*,*) 'DBG ZmatBuil - Treating Izm=4'
542 1 pfleura2
543 1 pfleura2
!     on a cree la premiere couche autour du premier centre
544 1 pfleura2
!     reste a finir les autres couches
545 1 pfleura2
      IndAFaire=1
546 1 pfleura2
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
547 1 pfleura2
      Do WHILE (Cafaire(IndAFaire).ne.0)
548 1 pfleura2
         IAt0=Cafaire(IndAfaire)
549 1 pfleura2
         if (debug) WRITE(*,*) "IndAFaire,IAt0,Lies_CF(IAt0,0)", &
550 1 pfleura2
              IndAFaire,IAt0,Lies_CF(IAt0,0)
551 1 pfleura2
         if (Lies_CF(IAt0,0).ge.1) THEN
552 1 pfleura2
!     IAt0 est l'atome pour lequel on construit la couche suivante
553 1 pfleura2
!     le premier atome lie est particulier car il definit l'orientation
554 1 pfleura2
!     de ce fragment
555 1 pfleura2
            IAti=Lies_CF(IAt0,1)
556 1 pfleura2
      WRITE(IOOUT,*) 'IAti:',IAti
557 1 pfleura2
            IAtd=IAt0
558 1 pfleura2
      WRITE(IOOUT,*) 'IAtd:',IAtd
559 1 pfleura2
            IAtv=Lies_CP(IAt0,1)
560 1 pfleura2
      WRITE(IOOUT,*) 'IAtv:',IAtv
561 1 pfleura2
            IIAtdh=1
562 1 pfleura2
            IAtdh=Lies_CF(IAtv,IIatdh)
563 1 pfleura2
            DO WHILE (Iatdh.eq.IAt0)
564 1 pfleura2
               IIatdh=IIatdh+1
565 1 pfleura2
               IAtdh=Lies_CF(IAtv,IIatdh)
566 1 pfleura2
            END DO
567 1 pfleura2
            IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
568 1 pfleura2
      WRITE(IOOUT,*) 'IAtdh:',IAtdh
569 1 pfleura2
            IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
570 1 pfleura2
               Izm=Izm+1
571 1 pfleura2
               if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
572 1 pfleura2
                    izm,IAti,IAtd,IAtv,IAtdh
573 1 pfleura2
               Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
574 1 pfleura2
                    ind_zmat,val_zmat ,x,y,z)
575 1 pfleura2
               DejaFait(IAti)=.TRUE.
576 1 pfleura2
            END IF
577 1 pfleura2
578 1 pfleura2
            CAfaire(IndFin)=Lies_CF(IAt0,1)
579 1 pfleura2
            IndFin=IndFin+1
580 1 pfleura2
!     Le traitement des autres est plus facile
581 1 pfleura2
            IAtdh=Lies_CF(IAt0,1)
582 1 pfleura2
            DO IAta=2,Lies_CF(IAt0,0)
583 1 pfleura2
               IAtI=Lies_CF(IAt0,IAta)
584 1 pfleura2
               CAfaire(IndFin)=IAtI
585 1 pfleura2
               IndFin=IndFin+1
586 1 pfleura2
587 1 pfleura2
               IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
588 1 pfleura2
                  Izm=Izm+1
589 1 pfleura2
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh,ind_zmat &
590 1 pfleura2
                       ,val_zmat,x,y,z)
591 1 pfleura2
                  DejaFait(IAti)=.TRUE.
592 1 pfleura2
               END IF
593 1 pfleura2
            END DO
594 1 pfleura2
            CAfaire(IndFin)=0
595 1 pfleura2
         END IF
596 1 pfleura2
         IndAFaire=IndAFaire+1
597 1 pfleura2
         if (debug) WRITE(IOOUT,*) "IndAfaire,CaFaire:",IndAfaire, &
598 1 pfleura2
            CaFaire
599 1 pfleura2
      END DO
600 1 pfleura2
601 1 pfleura2
602 1 pfleura2
!     On a fini de creer la molecule autour du premier atome 'centre'
603 1 pfleura2
!     Pour les autres c'est presque id... sauf que les axes sont deja fixes
604 1 pfleura2
!     On pourrait imaginer de mettre des atomes fictifs pour redefinir des axes
605 1 pfleura2
!     locaux... dans une autre version
606 1 pfleura2
!     V2.0
607 1 pfleura2
 9753 FFirst=.FALSE.
608 1 pfleura2
      if (debug) THEN
609 1 pfleura2
         WRITE(*,*) 'NbAt0r',NbAt0r,AtCP0
610 1 pfleura2
      END IF
611 1 pfleura2
      DO I=1,NbAt0r
612 1 pfleura2
!     Boucle pour trouver celui qui a le plus d'atomes CF
613 1 pfleura2
         IdAt0=0
614 1 pfleura2
         VCF=-1
615 1 pfleura2
         DO ICP0=1,NbAt0
616 1 pfleura2
            if (AtCP0(ICP0).NE.0) THEN
617 1 pfleura2
               ICat=AtCP0(ICP0)
618 1 pfleura2
               if (Lies_CF(ICAt,0).gt.VCF) THEN
619 1 pfleura2
                  IdAt0=ICat
620 1 pfleura2
                  IdxAt0=ICP0
621 1 pfleura2
                  VCF=Lies_CF(ICAt,0)
622 1 pfleura2
               END IF
623 1 pfleura2
            END IF
624 1 pfleura2
         END DO
625 1 pfleura2
         AtCP0(IdxAt0)=0
626 1 pfleura2
627 1 pfleura2
         IF (debug) WRITE(*,*) 'Adding fragment centered on ',IdAt0, &
628 1 pfleura2
              LIAISONS(IdAt0,0)
629 1 pfleura2
         IF (debug) WRITE(*,*) 'Center linked to ',Lies_CF(IdAt0,0), &
630 1 pfleura2
              ' atoms'
631 1 pfleura2
!     if (LIAISONS(IdAt0,0).EQ.0) goto 12345
632 1 pfleura2
633 1 pfleura2
         if (debug) THEN
634 1 pfleura2
            WRITE(IOOUT,*) 'LIAISONS tralala',nbcycle
635 1 pfleura2
            DO IAt=1,na
636 1 pfleura2
               WRITE(IOOUT,1003) iat,(LIAISONS(Iat,Jat),Jat=0,NMaxL)
637 1 pfleura2
            END DO
638 1 pfleura2
         END IF
639 1 pfleura2
 1003    FORMAT(1X,I5,' - ',12(I5))
640 1 pfleura2
!     Boucle pour voir quel est l'atome du fragment precedent qui est le plus
641 1 pfleura2
!     proche de celui ci
642 1 pfleura2
!     Rq: si ce sont des cycles lies, il pourrait etre malin de rechercher
643 1 pfleura2
!     a quoi il etait lie.
644 1 pfleura2
         norm2=1.e6
645 1 pfleura2
         Idx=0
646 1 pfleura2
         DO ICAt=1,Izm
647 1 pfleura2
            CALL vecteur(Idat0,ind_zmat(icat,1),x,y,z,vx1,vy1,vz1 &
648 1 pfleura2
                 ,norm1)
649 1 pfleura2
            if (norm2.ge.norm1) THEN
650 1 pfleura2
               norm2=norm1
651 1 pfleura2
               Idx=Icat
652 1 pfleura2
            END IF
653 1 pfleura2
         END DO
654 1 pfleura2
         IF (debug) WRITE(*,*) 'Link between fragments ',IdAt0, &
655 1 pfleura2
              ind_zmat(Idx,1), ' -- Idx=',Idx
656 1 pfleura2
657 1 pfleura2
! Despite the fact that this atom has officialy no CP atom
658 1 pfleura2
! we add one into its list; the one just found
659 1 pfleura2
         Lies_CP(IdAt0,0)=Lies_CP(IdAt0,0)+1
660 1 pfleura2
         Lies_CP(Idat0,Lies_CP(IdAt0,0))= ind_zmat(Idx,1)
661 1 pfleura2
         Lies_CP(Idat0,Lies_CP(IdAt0,0)+1)=0
662 1 pfleura2
663 1 pfleura2
664 1 pfleura2
!     on a l'indice... on va rajouter cet atome a la liste !
665 1 pfleura2
         izm=izm+1
666 1 pfleura2
         n1=ind_zmat(Idx,1)
667 1 pfleura2
!     Petit probleme si Idx<=2
668 1 pfleura2
         IF (Idx.EQ.1) THEN
669 1 pfleura2
!     Pb non negligeable ...
670 1 pfleura2
            IF (izm.ge.2) THEN
671 1 pfleura2
               IAtv=ind_zmat(2,1)
672 1 pfleura2
               IAtdh=ind_zmat(3,1)
673 1 pfleura2
            ELSE
674 1 pfleura2
               WRITE(*,*) '2 frag, le 1er a moins de 2 atomes..'
675 1 pfleura2
               WRITE(*,*) "J'ai merde... "
676 1 pfleura2
               STOP
677 1 pfleura2
            END IF
678 1 pfleura2
         ELSEIF (Idx.EQ.2) THEN
679 1 pfleura2
            IAtv=1
680 1 pfleura2
            IF (izm.ge.2) THEN
681 1 pfleura2
               IAtdh=ind_zmat(3,1)
682 1 pfleura2
            ELSE
683 1 pfleura2
               WRITE(*,*) '2 frag, le 1er a moins de 2 atomes...'
684 1 pfleura2
               WRITE(*,*) "J'ai merde... "
685 1 pfleura2
               STOP
686 1 pfleura2
            END IF
687 1 pfleura2
         ELSE
688 1 pfleura2
            IAtv=ind_zmat(Idx,2)
689 1 pfleura2
            IAtdh=ind_zmat(Idx,3)
690 1 pfleura2
         END IF
691 1 pfleura2
692 1 pfleura2
!     WRITE(*,*) "Outchilou IAtv, IAtdh", IAtv, IAtdh,n1
693 1 pfleura2
         CALL vecteur(n1,IdAt0,x,y,z,vx1,vy1,vz1,norm1)
694 1 pfleura2
         CALL vecteur(n1,IAtv,x,y,z,vx2,vy2,vz2,norm2)
695 1 pfleura2
         val=angle(vx1,vy1,vz1,norm1, &
696 1 pfleura2
              vx2,vy2,vz2,norm2)
697 1 pfleura2
         if (debug) WRITE(*,*) 'Angle val:1',IdAt0,n1,IAtv,val
698 1 pfleura2
         if ((abs(val).LE.10.).OR.(abs(pi-val).le.10.)) THEN
699 1 pfleura2
            IAtv=IAtdh
700 1 pfleura2
!     Ceci pose pb si Idx<=3... a traiter plus tard
701 1 pfleura2
            IF (Idx.ge.3) THEN
702 1 pfleura2
               IAtdh=ind_zmat(Idx,4)
703 1 pfleura2
            ELSE
704 1 pfleura2
               if (izm.ge.4) THEN
705 1 pfleura2
                  IAtdh=4
706 1 pfleura2
               ELSE
707 1 pfleura2
                  WRITE(*,*) 'Pb Angle Val+idx=3+izm<=4'
708 1 pfleura2
                  STOP
709 1 pfleura2
               END IF
710 1 pfleura2
            END IF
711 1 pfleura2
         END IF
712 1 pfleura2
         Call AddLigne(Izm,IdAt0,ind_zmat(Idx,1),IAtv,IAtdh, &
713 1 pfleura2
              ind_zmat,val_zmat,x,y,z)
714 1 pfleura2
         DejaFait(IdAt0)=.TRUE.
715 1 pfleura2
716 1 pfleura2
         IndFin=1
717 1 pfleura2
         IAtd=IdAt0
718 1 pfleura2
         n1=IAtdh
719 1 pfleura2
         IAtdh=IAtv
720 1 pfleura2
         IAtv=ind_zmat(Idx,1)
721 1 pfleura2
!     On ajoute les atomes lies a cet atome dans la liste a faire
722 1 pfleura2
         Do iaf=1,Lies_CF(IdAt0,0)
723 1 pfleura2
            IatI=Lies_CF(IdAt0,Iaf)
724 1 pfleura2
!     We check that this atom is not the one that is the closest to
725 1 pfleura2
!     the center atom
726 1 pfleura2
            if (IAtv.Ne.IAtI) THEN
727 1 pfleura2
               IF (debug) WRITE(*,*) 'Adding atom',IAtI, &
728 1 pfleura2
                    'to CAFaire in pos',iaf
729 1 pfleura2
               CAfaire(IndFin)=IAtI
730 1 pfleura2
               IndFin=IndFin+1
731 1 pfleura2
!     On les ajoute aussi dans la zmat
732 1 pfleura2
               If (ListAt(IAtI).AND.(.NOT.DejaFait(IatI))) THEN
733 1 pfleura2
                  izm=izm+1
734 1 pfleura2
       WRITE(*,*) "Outchili IAtv, IAtdh,IAtI", IAtv, IAtdh,IAtI
735 1 pfleura2
                  CALL vecteur(IAtI,IAtD,x,y,z,vx1,vy1,vz1,norm1)
736 1 pfleura2
                  CALL vecteur(IAtI,IAtv,x,y,z,vx2,vy2,vz2,norm2)
737 1 pfleura2
                  val=angle(vx1,vy1,vz1,norm1, vx2,vy2,vz2,norm2)
738 1 pfleura2
                  if (debug) &
739 1 pfleura2
                       WRITE(*,*) 'Angle val:2',IAtI,IAtD,IAtv,val
740 1 pfleura2
                  if ((abs(val).LE.10.).OR. &
741 1 pfleura2
                       (abs(180.-abs(val)).le.10.)) THEN
742 1 pfleura2
                     Call AddLigne(Izm,IAtI,IAtv,IAtdh,n1, &
743 1 pfleura2
                          ind_zmat,val_zmat,x,y,z)
744 1 pfleura2
                     DejaFait(Iati)=.TRUE.
745 1 pfleura2
                  ELSE
746 1 pfleura2
                     Call AddLigne(Izm,IAtI,IAtD,IAtv,IAtdh, &
747 1 pfleura2
                          ind_zmat,val_zmat,x,y,z)
748 1 pfleura2
                     DejaFait(Iati)=.TRUE.
749 1 pfleura2
                  END IF
750 1 pfleura2
               END IF
751 1 pfleura2
            END IF
752 1 pfleura2
         END DO
753 1 pfleura2
         CAfaire(IndFin)=0
754 1 pfleura2
755 1 pfleura2
!     On traite le reste de ce fragment !!
756 1 pfleura2
         IndAFaire=1
757 1 pfleura2
         WRITE(IOOUT,*) CaFaire
758 1 pfleura2
         Do WHILE (Cafaire(IndAFaire).ne.0)
759 1 pfleura2
            IAt0=Cafaire(IndAfaire)
760 1 pfleura2
            if (Lies_CF(IAt0,0).ge.1) THEN
761 1 pfleura2
!     IAt0 est l'atome pour lequel on construit la couche suivante
762 1 pfleura2
!     le premier atome lie est particulier car il definit l'orientation
763 1 pfleura2
!     de ce fragment
764 1 pfleura2
               Itmp=1
765 1 pfleura2
               IAti=Lies_CF(IAt0,Itmp)
766 1 pfleura2
               DO WHILE (DejaFait(IAti).AND. &
767 1 pfleura2
                    (Itmp.Le.Lies_CF(IAt0,0)))
768 1 pfleura2
                  Itmp=Itmp+1
769 1 pfleura2
                  IAti=Lies_CF(IAt0,ITmp)
770 1 pfleura2
               END DO
771 1 pfleura2
               If (ListAt(IAtI).AND.(.NOT.DejaFait(IatI))) THEN
772 1 pfleura2
                  IF (debug) WRITE(IOOUT,*) 'IAti:',IAti
773 1 pfleura2
                  IAtd=IAt0
774 1 pfleura2
                  IF (debug)     WRITE(IOOUT,*) 'IAtd:',IAtd
775 1 pfleura2
                  IAtv=Lies_CP(IAt0,1)
776 1 pfleura2
                  IF (debug)     WRITE(IOOUT,*) 'IAtv:',IAtv
777 1 pfleura2
                  IIAtdh=1
778 1 pfleura2
                  IAtdh=Lies_CF(IAtv,IIatdh)
779 1 pfleura2
                  DO WHILE (Iatdh.eq.IAt0)
780 1 pfleura2
                     IIatdh=IIatdh+1
781 1 pfleura2
                     IAtdh=Lies_CF(IAtv,IIatdh)
782 1 pfleura2
                  END DO
783 1 pfleura2
                  IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
784 1 pfleura2
                  IF (debug) WRITE(IOOUT,*) 'IAtdh:',IAtdh
785 1 pfleura2
                  Izm=Izm+1
786 1 pfleura2
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
787 1 pfleura2
                       ind_zmat,val_zmat,x,y,z)
788 1 pfleura2
                     DejaFait(Iati)=.TRUE.
789 1 pfleura2
!     Le traitement des autres est plus facile
790 1 pfleura2
                  IAtdh=Lies_CF(IAt0,ITmp)
791 1 pfleura2
                  DO IAta=ITmp+1,Lies_CF(IAt0,0)
792 1 pfleura2
                     IAtI=Lies_CF(IAt0,IAta)
793 1 pfleura2
                     If (ListAt(IAtI).AND.(.NOT.DejaFait(IatI))) &
794 1 pfleura2
                          THEN
795 1 pfleura2
                        Izm=Izm+1
796 1 pfleura2
                        Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
797 1 pfleura2
                             ind_zmat,val_zmat,x,y,z)
798 1 pfleura2
                        DejaFait(Iati)=.TRUE.
799 1 pfleura2
                     END IF
800 1 pfleura2
                  END DO
801 1 pfleura2
               END IF
802 1 pfleura2
803 1 pfleura2
!     On ajoute les atomes lies a cet atome dans la liste a faire
804 1 pfleura2
               Do iaf=1,Lies_CF(IAt0,0)
805 1 pfleura2
                  CAfaire(IndFin)=Lies_CF(IAt0,Iaf)
806 1 pfleura2
                  IndFin=IndFin+1
807 1 pfleura2
               END DO
808 1 pfleura2
               CAfaire(IndFin)=0
809 1 pfleura2
            END IF
810 1 pfleura2
            IndAFaire=IndAFaire+1
811 1 pfleura2
         END DO
812 1 pfleura2
12345    CONTINUE
813 1 pfleura2
      END DO
814 1 pfleura2
815 1 pfleura2
      WRITE(*,*) 'TOTOTOTOTOTOTOTOT'
816 1 pfleura2
817 1 pfleura2
      if (debug) THEN
818 1 pfleura2
         WRITE(*,*) 'ind_zmat'
819 1 pfleura2
         DO I=1,na
820 1 pfleura2
            WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
821 1 pfleura2
         END DO
822 1 pfleura2
      END IF
823 1 pfleura2
824 1 pfleura2
      IndFin=1
825 1 pfleura2
      DO I=1,Izm0
826 1 pfleura2
         Iat=ind_zmat(I,1)
827 1 pfleura2
         Do iaf=1,Lies_CF(Iat,0)
828 1 pfleura2
            IIat=Lies_CF(Iat,iaf)
829 1 pfleura2
            if (ListAt(iIAt).AND. &
830 1 pfleura2
                 (.NOT.DejaFait(iIAt)))  THEN
831 1 pfleura2
               IAti=IIat
832 1 pfleura2
               WRITE(IOOUT,*) 'IAti:',IAti
833 1 pfleura2
               IAtd=IAt0
834 1 pfleura2
               WRITE(IOOUT,*) 'IAtd:',IAtd
835 1 pfleura2
               IAtv=Lies_CP(IAt0,1)
836 1 pfleura2
               WRITE(IOOUT,*) 'IAtv:',IAtv
837 1 pfleura2
               IIAtdh=1
838 1 pfleura2
               IAtdh=Lies_CF(IAtv,IIatdh)
839 1 pfleura2
               DO WHILE (Iatdh.eq.IAt0)
840 1 pfleura2
                  IIatdh=IIatdh+1
841 1 pfleura2
                  IAtdh=Lies_CF(IAtv,IIatdh)
842 1 pfleura2
               END DO
843 1 pfleura2
               IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
844 1 pfleura2
               WRITE(IOOUT,*) 'IAtdh:',IAtdh
845 1 pfleura2
               IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
846 1 pfleura2
                  Izm=Izm+1
847 1 pfleura2
                  if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
848 1 pfleura2
                       izm,IAti,IAtd,IAtv,IAtdh
849 1 pfleura2
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
850 1 pfleura2
                       ind_zmat,val_zmat ,x,y,z)
851 1 pfleura2
                  DejaFait(IAti)=.TRUE.
852 1 pfleura2
                  CaFaire(IndFin)=iIat
853 1 pfleura2
                  IndFin=IndFin+1
854 1 pfleura2
               END IF
855 1 pfleura2
            END IF
856 1 pfleura2
         END DO
857 1 pfleura2
      END DO
858 1 pfleura2
      CAfaire(IndFin)=0
859 1 pfleura2
860 1 pfleura2
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
861 1 pfleura2
862 1 pfleura2
      if (debug) WRITE(*,*) 'DBG ZmatBuil - Toto'
863 1 pfleura2
      if (debug) THEN
864 1 pfleura2
         WRITE(*,*) 'ind_zmat'
865 1 pfleura2
         DO I=1,na
866 1 pfleura2
            WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
867 1 pfleura2
         END DO
868 1 pfleura2
      END IF
869 1 pfleura2
870 1 pfleura2
!     on a cree la premiere couche autour du premier centre
871 1 pfleura2
!     reste a finir les autres couches
872 1 pfleura2
      IndAFaire=1
873 1 pfleura2
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
874 1 pfleura2
      Do WHILE (Cafaire(IndAFaire).ne.0)
875 1 pfleura2
         IAt0=Cafaire(IndAfaire)
876 1 pfleura2
         if (debug) WRITE(*,*) "IndAFaire,IAt0,Lies_CF(IAt0,0)", &
877 1 pfleura2
              IndAFaire,IAt0,Lies_CF(IAt0,0)
878 1 pfleura2
         if (Lies_CF(IAt0,0).ge.1) THEN
879 1 pfleura2
!     IAt0 est l'atome pour lequel on construit la couche suivante
880 1 pfleura2
!     le premier atome lie est particulier car il definit l'orientation
881 1 pfleura2
!     de ce fragment
882 1 pfleura2
            IAti=Lies_CF(IAt0,1)
883 1 pfleura2
      WRITE(IOOUT,*) 'IAti:',IAti
884 1 pfleura2
            IAtd=IAt0
885 1 pfleura2
      WRITE(IOOUT,*) 'IAtd:',IAtd
886 1 pfleura2
            IAtv=Lies_CP(IAt0,1)
887 1 pfleura2
      WRITE(IOOUT,*) 'IAtv:',IAtv
888 1 pfleura2
            IIAtdh=1
889 1 pfleura2
            IAtdh=Lies_CF(IAtv,IIatdh)
890 1 pfleura2
            DO WHILE (Iatdh.eq.IAt0)
891 1 pfleura2
               IIatdh=IIatdh+1
892 1 pfleura2
               IAtdh=Lies_CF(IAtv,IIatdh)
893 1 pfleura2
            END DO
894 1 pfleura2
            IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
895 1 pfleura2
      WRITE(IOOUT,*) 'IAtdh:',IAtdh
896 1 pfleura2
            IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
897 1 pfleura2
               Izm=Izm+1
898 1 pfleura2
               if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
899 1 pfleura2
                    izm,IAti,IAtd,IAtv,IAtdh
900 1 pfleura2
               Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
901 1 pfleura2
                    ind_zmat,val_zmat ,x,y,z)
902 1 pfleura2
               DejaFait(IAti)=.TRUE.
903 1 pfleura2
            END IF
904 1 pfleura2
905 1 pfleura2
            CAfaire(IndFin)=Lies_CF(IAt0,1)
906 1 pfleura2
            IndFin=IndFin+1
907 1 pfleura2
!     Le traitement des autres est plus facile
908 1 pfleura2
            IAtdh=Lies_CF(IAt0,1)
909 1 pfleura2
            DO IAta=2,Lies_CF(IAt0,0)
910 1 pfleura2
               IAtI=Lies_CF(IAt0,IAta)
911 1 pfleura2
               CAfaire(IndFin)=IAtI
912 1 pfleura2
               IndFin=IndFin+1
913 1 pfleura2
914 1 pfleura2
               IF (ListAt(IAti).AND.(.NOT.DejaFait(IAti))) THEN
915 1 pfleura2
                  IF (debug) &
916 1 pfleura2
         WRITE(*,*) 'DBG ZmatBuil Toto adding',Iati,Izm
917 1 pfleura2
                  Izm=Izm+1
918 1 pfleura2
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh,ind_zmat &
919 1 pfleura2
                       ,val_zmat,x,y,z)
920 1 pfleura2
                  DejaFait(IAti)=.TRUE.
921 1 pfleura2
               END IF
922 1 pfleura2
            END DO
923 1 pfleura2
            CAfaire(IndFin)=0
924 1 pfleura2
         END IF
925 1 pfleura2
         IndAFaire=IndAFaire+1
926 1 pfleura2
         if (debug) WRITE(IOOUT,*) "Toto IndAfaire,CaFaire:",IndAfaire, &
927 1 pfleura2
            CaFaire
928 1 pfleura2
      END DO
929 1 pfleura2
930 1 pfleura2
931 1 pfleura2
      if (debug) THEN
932 1 pfleura2
         WRITE(*,*) 'ind_zmat'
933 1 pfleura2
         DO I=1,na
934 1 pfleura2
            WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
935 1 pfleura2
         END DO
936 1 pfleura2
      END IF
937 1 pfleura2
938 1 pfleura2
      END