Statistiques
| Révision :

root / src / Calc_zmat.f90 @ 8

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

1
      SUBROUTINE Calc_Zmat(na,atome,ind_zmat,val_zmat,x,y,z, &
2
           r_cov,  fact)
3

    
4
      use Path_module, only : Pi,max_Z, NMaxL, Nom, KINT, KREAL
5
      use Io_module, only :  IoOut
6

    
7
      IMPLICIT NONE
8

    
9
! Number of atoms
10
      integer(KINT), INTENT(IN) :: na
11
! Mass number of atoms
12
      INTEGER(KINT), INTENT(IN) :: atome(na)
13
! Coordinates
14
      real(KREAL), INTENT(IN) ::  x(na),y(na),z(na)
15
! Covalent radii
16
      REAL(KREAL), INTENT(IN) :: R_cov(Max_Z)
17
! Factor to determine connectivity
18
      REAL(KREAL) :: Fact
19
! Zmat index and values
20
      integer(KINT), INTENT(OUT) :: ind_zmat(na,5)
21
      real(KREAL),INTENT(OUT) ::  val_zmat(na,3)
22

    
23

    
24
      INTEGER(KINT), ALLOCATABLE :: LIAISONS(:,:) !(na,0:NMaxL
25
      INTEGER(KINT), ALLOCATABLE :: LiaisonsBis(:,:) ! na,0:NMaxL
26
      INTEGER(KINT), ALLOCATABLE :: LieS_CP(:,:),LieS_CF(:,:) ! (na,0:NMaxL)
27
      INTEGER(KINT), ALLOCATABLE :: NbAt3(:,:) !(Na,2)
28
      INTEGER(KINT), ALLOCATABLE :: CycleAt(:) !(Na)
29
      INTEGER(KINT), ALLOCATABLE :: CAFaire(:) ! Na
30

    
31
      Integer(KINT) :: AtTypCycl(Max_Z)
32
      Integer(KINT) :: NbCycle
33
      real(KREAL) ::  vx1,vy1,vz1,norm1
34
      real(KREAL) ::  vx2,vy2,vz2,norm2
35
      real(KREAL) ::  val
36
      Logical PasFini,AtTerm,Debug,Bicycle,FLie3,FirstCycle
37
      LOGICAL, ALLOCATABLE ::  Former3(:), DejaFait(:) ! Na
38
      Logical FTmp
39
      LOGICAL F1213, F1223,F1323
40

    
41
      INTEGER(KINT) :: I,J,K, IZM, N1,n2, N3, IL, JL
42
      INTEGER(KINT) :: IdAt0, I3, IAt3, Idx3, I1, I0, IAf
43
      INTEGER(KINT) :: Izm1,Izm2, IZm3,IZm4, VCF, ICat
44
      INTEGER(KINT) :: IOld, IndFin, IdAtL, IndAFaire
45
      INTEGER(KINT) :: IAti, IAtD, IAtv, IIAtDh, IAtDh
46
      INTEGER(KINT) :: IAt0, IAta, IAt, Jat,Idx
47
      INTEGER(KINT) :: ITmp, IAtTmp, NbAt0
48
      INTEGER(KINT) :: NbLies, ICycle, IMax
49
      REAL(KREAL) :: d,d12, d13,d32
50

    
51
  INTERFACE
52
     function valid(string) result (isValid)
53
       CHARACTER(*), intent(in) :: string
54
       logical                  :: isValid
55
     END function VALID
56

    
57
     FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
58
        use Path_module, only :  Pi,KINT, KREAL
59
       real(KREAL) ::  v1x,v1y,v1z,norm1
60
       real(KREAL) ::  v2x,v2y,v2z,norm2
61
       real(KREAL) ::  angle
62
     END FUNCTION ANGLE
63

    
64

    
65
  END INTERFACE
66

    
67

    
68

    
69
      debug=valid("Calc_zmat")
70
      if (debug) WRITE(*,*) "========== Entering Calc_zmat =========="
71

    
72
      FirstCycle=.TRUE.
73
      FTmp=.FALSE.
74
      NbCycle=0
75
      DO i=1,na
76
         DO J=1,5
77
            Ind_Zmat(i,J)=0
78
         END DO
79
      END DO
80
      ALLOCATE(Former3(Na), DejaFait(Na))
81
      ALLOCATE(CaFaire(Na), CycleAt(Na))
82
      ALLOCATE(NbAt3(Na,2))
83
      ALLOCATE(Liaisons(na,0:NMaxL),LiaisonsBis(na,0:NMaxL))
84
      ALLOCATE(Lies_CP(na,0:NMaxL), Lies_CF(na,0:NMaxL))
85

    
86

    
87
      WRITE(IOOUT,*) Na
88
      WRITE(IOOUT,*) 'Calc_zmat'
89
      DO I=1,na
90
         WRITE(IOOUT,*) NOM(Atome(I)),X(i),y(i),z(i)
91
      END DO
92

    
93
      if (na.le.2) THEN
94
         WRITE(*,*) "I do not work for less than 2 atoms :-p"
95
         RETURN
96
      END IF
97

    
98
!     Cas particulier: 3 atomes ou moins...
99
      if (Na.eq.3) THEN
100
         d12=sqrt((x(1)-x(2))**2+(y(1)-y(2))**2+(z(1)-z(2))**2)
101
         d13=sqrt((x(1)-x(3))**2+(y(1)-y(3))**2+(z(1)-z(3))**2)
102
         d32=sqrt((x(3)-x(2))**2+(y(3)-y(2))**2+(z(3)-z(2))**2)
103
         F1213=(d12<=d13)
104
         F1323=(d13<=d32)
105
         F1223=(d12<=d32)
106
         if (debug) THEN
107
            WRITE(*,*) "DEBUG Calc_Zmat 3 atoms"
108
            WRITE(*,*) "d12,d13,d32:",d12,d13,d32
109
            WRITE(*,*) "F1213,F1323,F1223:",F1213,F1323,F1223
110
         END IF
111
         if (F1213) THEN
112
            if (F1323) THEN
113
               ind_zmat(1,1)=3
114
               ind_zmat(2,1)=1
115
               ind_zmat(2,2)=3
116
               ind_zmat(3,1)=2
117
               ind_zmat(3,2)=1
118
               ind_zmat(3,3)=3
119
            ELSE
120
               ind_zmat(1,1)=1
121
               ind_zmat(2,1)=2
122
               ind_zmat(2,2)=1
123
               ind_zmat(3,1)=3
124
               ind_zmat(3,2)=2
125
               ind_zmat(3,3)=1
126
            END IF
127
         ELSE
128
            IF (F1223) THEN
129
               ind_zmat(1,1)=2
130
               ind_zmat(2,1)=1
131
               ind_zmat(2,2)=2
132
               ind_zmat(3,1)=3
133
               ind_zmat(3,2)=1
134
               ind_zmat(3,3)=2
135
            ELSE
136
               ind_zmat(1,1)=2
137
               ind_zmat(2,1)=3
138
               ind_zmat(2,2)=2
139
               ind_zmat(3,1)=1
140
               ind_zmat(3,2)=3
141
               ind_zmat(3,3)=2
142
            END IF
143
         END IF
144
         IF (debug) THEN
145
            DO i=1,na
146
               WRITE(*,*) (ind_zmat(i,j),j=1,5)
147
            END DO
148
         END IF
149

    
150
!     We have ind_zmat, we fill val_zmat
151
         val_zmat(1,1)=0.
152
         val_zmat(1,2)=0.
153
         val_zmat(1,3)=0.
154
         n2=ind_zmat(2,1)
155
         n1=ind_zmat(2,2)
156
         d=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2+ &
157
              (z(n1)-z(n2))**2)
158
         val_zmat(2,1)=d
159
         val_zmat(2,2)=0.
160
         val_zmat(2,3)=0.
161
         n1=ind_zmat(3,1)
162
         n2=ind_zmat(3,2)
163
         n3=ind_zmat(3,3)
164
         CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
165
         if (debug) write(*,*) n1,n2,norm1
166
         val_zmat(3,1)=norm1
167

    
168
         CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
169
         val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
170
         if (debug) write(*,*) n2,n3,norm2,val
171
         val_zmat(3,2)=val
172
         val_zmat(3,3)=0.
173

    
174
         RETURN
175
      END IF
176

    
177

    
178
!     Premiere etape : calcul des connectivites
179
      DO I=1,na
180
         DejaFait(I)=.FALSE.
181
         Former3(I)=.False.
182
         Do J=0,NMaxL
183
            Liaisons(I,j)=0
184
            LiaisonsBis(I,j)=0
185
         END DO
186
         DO J=1,5
187
            ind_zmat(I,J)=0
188
         END DO
189
      END DO
190

    
191
      if (debug) WRITE(IOOUT,*) "Liaisons initialiazed"
192
!     DO IL=1,na
193
!     WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
194
!     WRITE(IOOUT,1003) Il,(LIAISONSBis(IL,JL),JL=0,NMaxL)
195
!     END DO
196

    
197
      Call CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)
198

    
199
      if (debug) THEN
200
         WRITE(IOOUT,*) "Connections calculated"
201
         DO IL=1,na
202
            WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
203
         END DO
204
      END IF
205

    
206

    
207
!     on va maintenant trier ces connectivites en 2 :
208
!     Lies_CF : vers l'exterieur de la molecule
209
!     Lies_CP : vers l'interieur de la molecule
210
!     Pour cela, on procede 'a l'envers' : on regarde les atomes
211
!     auxquels sont lies les atomes attaches -> on remplit Lies_CF/P
212
!     et on supprime ces atomes...
213

    
214
!     v2.0 on copie Liaison dans LiaisonBis. On analyse LiaisonBis
215
!     tout en modifiant Liaison. Cela evite de fausser l'analyse: un atome
216
!     qui est lie initialement a 2 atomes (et qui n'est donc pas terminal)
217
!     peut devenir terminal en milieu de traitement si on annule la liaison
218
!     qui le liait a un atome terminal situ? avant lui...
219
!     ex: H2O code dans l'ordre H,O,H.
220
      PasFini=.True.
221
      AtTerm=.True.
222
      DO WHILE (PasFini.AND.AtTerm)
223
         PasFini=.False.
224
         AtTerm=.False.
225
         DO I=1,na
226
            DO J=0,NMaxL
227
               LiaisonsBis(I,J)=Liaisons(I,J)
228
            END DO
229
         END DO
230
!     DO IL=1,na
231
!     WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
232
!     WRITE(IOOUT,1003) Il,(LIAISONSBis(IL,JL),JL=0,NMaxL)
233
!     END DO
234
!     WRITE(IOOUT,*) "=================================="
235

    
236

    
237
         DO I=1,na
238
            if ((LiaisonsBis(I,0).eq.1).AND.Liaisons(I,0).ne.0) THEN
239
               AtTerm=.True.
240
!     On enleve cet atome
241
               Liaisons(I,0)=0
242
               NbLies=Lies_CP(I,0)+1
243
               Lies_CP(I,NbLies)=Liaisons(I,1)
244
               Lies_CP(I,0)=NbLies
245
               Liaisons(Liaisons(I,1),0)=Liaisons(Liaisons(I,1),0)-1
246

    
247
               NbLies=Lies_CF(Liaisons(I,1),0)+1
248
               Lies_CF(Liaisons(I,1),NbLies)=I
249
               Lies_CF(Liaisons(I,1),0)=NbLies
250

    
251
               Call Annul(Liaisons,Liaisons(I,1),I)
252

    
253
            END IF
254

    
255
            if (Liaisons(I,0).ge.1) THEN
256
               PasFini=.TRUE.
257
            END IF
258

    
259
         END DO
260
         if (debug) WRITE(*,*) 'PasFini, AtTerm',PasFini,AtTerm
261
         If (PasFini.AND.(.Not.AtTerm)) THEN
262
!     Pas d'atomes terminaux lors de l'exploration precedente.
263
!     On a soit une erreur, soit un cycle. Je ne sais pas encore gerer
264
!     tous les  cas :  on affiche un warning
265
            WRITE(*,*) "Je ne trouve pas d'atomes terminaux"
266
            WRITE(*,*) "Possibilite de molecule cyclique"
267
            WRITE(*,*) "Cas en cours de traitement: verifier l'output"
268
!     On va d'abord voir si on a des atomes li?s a plus de 2 centres.
269
            AtTerm=.TRUE.
270
            Bicycle=.False.
271
            ICycle=0
272
            IMax=0
273
            DO I=1,na
274
               if (Liaisons(I,0).gt.2) THEN
275
                  ICycle=ICycle+1
276
                  BiCycle=.True.
277
                  If (Liaisons(I,0).gt.IMax) Imax=Liaisons(I,0)
278
                  Former3(I)=.True.
279
               END IF
280
            END DO
281
            IF (BiCycle) THEN
282
!     On a des atomes lies a 3 ou plus d'atomes...
283
!     donc soit des bicycles, soit plusieurs
284
!     cycles attaches par un sommet...
285
               WRITE(*,*) "On dirait un bicyle... on essaie"
286
               WRITE(*,*) ICycle, "atomes lies a plus de 2 atomes"
287
               WRITe(*,*) "Nb Attaches max:",IMax
288
!     Plusieurs cas possibles pour un bicycle:
289
!     1) 2 cycles relies par un sommet, ICycle=2, IMax=3
290
!     2) 2 cycles relies par une arrete, ICycle=2, Imax=3
291
               If (Imax.Eq.3) THEN
292
!     on doit pouvoir traiter celui la...
293
!     On classe les atomes li?s a trois atomes, par le nombre d'atomes
294
!     lies a trois atomes
295
!     auxquels ils sont li?s ... interessant pour les composes asterisques.
296
                  I3=0
297
                  FLie3=.False.
298
                  DO I=1,Na
299
                     If (Liaisons(I,0).EQ.3) THEN
300
                        I3=I3+1
301
                        K=0
302
                        DO J=1,Liaisons(I,0)
303
                           If (Liaisons(Liaisons(I,J),0).Eq.3) THEN
304
                              k=k+1
305
                              FLie3=.True.
306
                           END IF
307
                        END DO
308
                        NbAt3(I3,2)=k
309
                        NbAt3(I3,1)=I
310
                     END IF
311
                  END DO
312
!     A-t-on deux atomes a 3 lies ensemble ?
313
                  IF (FLie3) THEN
314
!     on les classe pas nb at lies a 3
315
                     IAt3=0
316
                     Idx3=0
317
                     DO I=1,I3
318
                        IF (NbAt3(I,2).ge.Iat3) THEN
319
                           IAt3=NbAt3(I,2)
320
                           Idx3=NbAt3(I,1)
321
                        END IF
322
                     END DO
323
!     On va enlever ces liaisons entre atomes a 3, en les mettant
324
!     en CF de l'atome 'central'
325
                     if (debug) WRITE(*,*) "Atome ",Idx3, &
326
                          " retenu, li? a",IAt3," atomes 3"
327
                     DO I=1,Liaisons(Idx3,0)
328
                        I1=Liaisons(Idx3,I)
329
                        If (Liaisons(I1,0).EQ.3) THEN
330
                           if (debug) WRITE(*,*) "Traitement ",I1,Idx3
331
                           NbLies=Lies_CF(Idx3,0)+1
332
                           Lies_CF(Idx3,NbLies)=I1
333
                           Lies_CF(Idx3,0)=NbLies
334
                           NbLies=Lies_CP(I1,0)+1
335
                           Lies_CP(I1,NbLies)=Idx3
336
                           Call Annul(Liaisons,I1,Idx3)
337
                           Call Annul(Liaisons,Idx3,I1)
338
                           Liaisons(Idx3,0)=Liaisons(Idx3,0)-1
339
                           Liaisons(I1,0)=Liaisons(I1,0)-1
340
                        END IF
341
                     END DO
342
                  ELSE
343
                     WRITE(*,*) "Aucune liaisons entre deux atomes li?s"
344
                     WRITE(*,*) "Pas encore trait?"
345
                     STOP
346
                  END IF        !FLie3=T ?
347
               END IF           !Imax=3 ?
348
            ELSE                ! BiCyle ?
349
!     Un seul cycle. Facile :-)
350
               WRITE(*,*) "Un  cycle identifie..."
351
               NbCycle=NbCycle+1
352
               if (debug) WRITE(*,*) 'Considering cycle ',NbCycle
353
               If (NbCycle.ge.2) THEN
354
!     si ce n'est pas le premier cycle que l'on met, on regarde si parmi les
355
!     atomes du cycle, l'un d'eux etait avant attache a 3 atomes...
356
                  I=1
357
                  DO WHILE (Liaisons(I,0).ne.2)
358
                     I=I+1
359
                  END DO
360
                  if (debug) WRITE(*,*) "I>2:",I,Former3(I)
361
                  FTmp=Former3(I)
362
                  I0=I
363
                  IOld=I
364
                  I1=I
365
                  I=Liaisons(I,1)
366
                  DO WHILE (.NOT.FTmp)
367
                     I1=Liaisons(I,1)
368
                     If (I1.Eq.IOld) I1=Liaisons(I,2)
369
                     FTmp=Former3(I1)
370
                     if (debug) WRITE(*,*) "I,I1,I0,IOld",I,I1,I0, &
371
                          IOld,FTmp
372
                     IF (I1.eq.I0) FTmp=.TRUE.
373
                     if (debug) WRITE(*,*) "I,I1,I0,IOld",I,I1,I0, &
374
                          IOld,FTmp
375
                     IOld=I
376
                     I=I1
377
                  END DO
378
                  if (debug) WRITE(*,*)  "I,I1,I0,IOld",I,I1,I0, &
379
                       IOld,FTmp
380
                  IF (Former3(I1)) THEN
381
!     On a notre atome particulier ... cool :-)
382
                     if (debug) WRITE(*,*) "Atome former3",I1
383
                     ITmp=I1
384
                  END IF
385
               END IF           ! NbCycle >=2
386
               IF (.NOT.Ftmp) THEN
387
                  if (debug) WRITE(*,*) "Pas trouve de former3"
388
!     on regarde si il y a un atome particulier.. ie un heteroatome ou autre.
389
                  DO I=1,Max_Z
390
                     AtTypCycl(I)=0
391
                  END DO
392
                  DO I=1,na
393
                     if (Liaisons(I,0).eq.2) &
394
                          AtTypCycl(Atome(I))=AtTypCycl(Atome(I))+1
395
                     if (debug) WRITE(*,*) I,Atome(I), &
396
                          AtTypCycl(Atome(I)), &
397
                          Liaisons(I,0)
398
                  END DO
399
                  Itmp=NmaxL+1
400
                  IAtTmp=0
401
                  DO I=1,Max_Z
402
                     If ((AtTypCycl(I).gt.0).and.(AtTypCycl(I).le.Itmp)) &
403
                          THEN
404
                        Itmp=AtTypCycl(I)
405
                        IAtTmp=I
406
                     END IF
407
                  END DO
408
                  if (debug) WRITE(*,*) 'Itmp,IAtTmp:',Itmp,IAtTmp
409
!     On a le type de l'atome particulier... on va prendre le premier venu
410
                  DO I=na,1,-1
411
                     If ((Atome(I).eq.IAtTmp).AND.(Liaisons(I,0).Eq.2)) &
412
                          Itmp=I
413
                  END DO
414
               END IF
415
               if (debug) WRITE(*,*) 'Atome ',Itmp,'(',IAtTmp,') retenu'
416
!     On va tricher en enlevant les liaisons de cet atome a la main...
417
               I0=Itmp
418
               I1=Liaisons(I0,1)
419
               DO WHILE (I1.Ne.ITmp)
420
                  if (debug) WRITE(*,*) "Going from",I0,"to ",I1
421
!     On rajoute ces deux liaisons en CF de Itmp
422
                  NbLies=Lies_CF(I0,0)+1
423
                  Lies_CF(I0,NbLies)=Liaisons(I0,1)
424
                  Lies_CF(I0,0)=NbLies
425
!     et les liaisons vers Itmp deviennent des CP pour les autres atomes
426
                  NbLies=Lies_CP(I1,0)+1
427
                  Lies_CP(I1,NbLies)=I0
428
                  Lies_CP(I1,0)=NbLies
429

    
430
                  Call Annul(Liaisons,I1,I0)
431
                  Liaisons(I1,0)=Liaisons(I1,0)-1
432
                  Liaisons(I0,0)=Liaisons(I0,0)-1
433

    
434
                  DO IL=1,na
435
                     WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
436
                  END DO
437
                  I0=I1
438
                  I1=Liaisons(I0,1)
439

    
440
               end do
441
               Call Annul(Liaisons,I1,I0)
442
               Liaisons(I1,0)=Liaisons(I1,0)-1
443
               Liaisons(I0,0)=Liaisons(I0,0)-1
444
            END IF
445
         END IF
446
         DO IL=1,na
447
            WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
448
         END DO
449
!     WRITE(IOOUT,*) "=================================="
450
!     WRITE(IOOUT,*) "=================================="
451
      END DO
452

    
453
!     WRITE(IOOUT,*) "=================================="
454
!     Il ne reste plus que des atomes lies a rien...
455
!     ce sont les 'centres' de la molecule. On repart
456
!     d'eux pour construire le reste de la molecule.
457

    
458
 1003 FORMAT(1X,I4,' - ',15(I5))
459
!     Avant de partir, on va classer, pour chaque atome, les atomes CF par
460
!     nombre de masse decroissant.
461
      DO I=1,na
462
!     WRITE(*,*) I,Lies_CF(I,0),(Atome(Lies_CF(I,J)),j=1,Lies_CF(I,0))
463
         DO J=1,Lies_CF(I,0)-1
464
            DO K=J+1,Lies_CF(I,0)
465
               if (atome(Lies_CF(I,K))>atome(Lies_CF(I,J))) THEN
466
                  Itmp=Lies_CF(I,K)
467
                  Lies_CF(I,K)=Lies_CF(I,J)
468
                  Lies_CF(I,J)=Itmp
469
               END IF
470
            END DO
471
         END DO
472
!     WRITE(*,*) I,Lies_CF(I,0),(Atome(Lies_CF(I,J)),j=1,Lies_CF(I,0))
473
      END DO
474

    
475
      IF (Debug) THEN
476
         WRITE(IOOUT,*) 'LIAISONS'
477
         DO I=1,na
478
            WRITE(IOOUT,1003) i,(LIAISONS(I,J),J=0,NMaxL)
479
         END DO
480

    
481
         WRITE(IOOUT,*) 'LIes_CF'
482
         DO I=1,na
483
            WRITE(IOOUT,1003) i,(LIeS_CF(I,J),J=0,NMaxL)
484
         END DO
485

    
486
         WRITE(IOOUT,*) 'LIes_CP'
487
         DO I=1,na
488
            WRITE(IOOUT,1003) i,(LIeS_CP(I,J),J=0,NMaxL)
489
         END DO
490
      END IF
491

    
492

    
493
!     On compte le nb d'atomes sans atomes CP (centripetes)
494
      NbAt0=0
495
      DO I=1,na
496
         IF (Lies_CP(I,0).eq.0) THEN
497
            NbAt0=NbAt0+1
498
            IF (DEBUG) WRITE(*,*) "Center atom",NbAt0,I,Atome(I)
499
         END IF
500
      END DO
501

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

    
506
      IZm=1
507
!     Boucle pour trouver celui qui a le plus d'atomes CF
508
      IdAt0=0
509
      VCF=0
510
      DO ICAt=1,na
511
         if ((Lies_CP(ICaT,0).eq.0).AND.(Lies_CF(ICAt,0).gt.VCF)) THEN
512
            IdAt0=ICat
513
            VCF=Lies_CF(ICAt,0)
514
         END IF
515
      END DO
516
      Lies_CP(IdAt0,0)=-1
517

    
518
!     WRITE(IOOUT,*) 'Atome 1', IdAt0, Nom(Atome(IdAt0))
519
      Izm1=IdAt0
520
      ind_zmat(izm,1)=Izm1
521
      ind_zmat(izm,2)=0
522
      ind_zmat(izm,3)=0
523
      ind_zmat(izm,4)=0
524
      ind_zmat(izm,5)=0
525
      val_zmat(izm,1)=0.0
526
      val_zmat(izm,2)=0.0
527
      val_zmat(izm,3)=0.0
528
      DejaFait(Izm1)=.True.
529

    
530
!     Les atomes CF lies a IdAt0 sont mis en attente pour
531
!     etre traites
532

    
533
      IndFin=1
534
      Do iaf=1,Lies_CF(IdAt0,0)
535
         CAfaire(IndFin)=Lies_CF(IdAt0,Iaf)
536
         IndFin=IndFin+1
537
      END DO
538
      CAfaire(IndFin)=0
539

    
540
!     On construit la premiere couche qui est speciale car elle fixe les
541
!     axes.
542
!     Plusieurs cas sont possibles suivant le nb d'atomes CF
543

    
544
      IF (Lies_CF(IdAt0,0).ge.3) THEN
545
         if (debug) WRITE(*,*) 'Cas tres simple,',IdAt0, &
546
              ' li? a plus de 3 atomes'
547
         IZm2=Lies_CF(IdAt0,1)
548
         IZm3=Lies_CF(IdAt0,2)
549
!     WRITE(IOOUT,*) 'Atome 2' , Izm2, Nom(Atome(Izm2))
550
!     WRITE(IOOUT,*) 'Atome 3' ,Izm3, Nom(Atome(izm3))
551

    
552
         Izm=2
553
         CALL vecteur(IZm1,IZm2,x,y,z,vx1,vy1,vz1,norm1)
554
!     write(*,11) n1,n2,norm1
555

    
556
         ind_zmat(izm,1)=izm2
557
         ind_zmat(izm,2)=izm1
558
         ind_zmat(izm,3)=0
559
         ind_zmat(izm,4)=0
560
         ind_zmat(izm,5)=0
561
         val_zmat(izm,1)=norm1
562
         val_zmat(izm,2)=0.0
563
         val_zmat(izm,3)=0.0
564
         DejaFait(Izm2)=.TRUE.
565

    
566
         Izm=3
567
         CALL vecteur(izm1,izm3,x,y,z,vx2,vy2,vz2,norm2)
568
         val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
569

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

    
572
         ind_zmat(izm,1)=izm3
573
         ind_zmat(izm,2)=izm1
574
         ind_zmat(izm,3)=izm2
575
         ind_zmat(izm,4)=0
576
         ind_zmat(izm,5)=0
577
         val_zmat(izm,1)=norm2
578
         val_zmat(izm,2)=val
579
         val_zmat(izm,3)=0.0
580
         DejaFait(Izm3)=.TRUE.
581

    
582
         DO IdAtl=3,Lies_CF(IdAt0,0)
583
            Izm4= Lies_CF(IdAt0,IdAtl)
584
!     WRITE(IOOUT,*) 'Atome ',Izm+1,Izm4,Nom(Atome(Izm4)),Izm1,
585
!     $Izm2,Izm3
586
            Izm=Izm+1
587
            Call AddLigne(Izm,Izm4,Izm1,Izm2,Izm3,ind_zmat,val_zmat, &
588
                 x,y,z)
589
            DejaFait(Izm4)=.TRUE.
590
         END DO
591
      END IF
592

    
593
      IF (Lies_CF(Izm1,0).eq.2) THEN
594
         if (debug) WRITE(*,*) 'Cas simple,',IdAt0, &
595
              ' li? a 2 atomes'
596

    
597
         IZm2=Lies_CF(IdAt0,1)
598
         IZm3=Lies_CF(IdAt0,2)
599
         Izm=2
600
         CALL vecteur(IZm1,IZm2,x,y,z,vx1,vy1,vz1,norm1)
601
!     write(*,11) n1,n2,norm1
602

    
603
         ind_zmat(izm,1)=izm2
604
         ind_zmat(izm,2)=izm1
605
         ind_zmat(izm,3)=0
606
         ind_zmat(izm,4)=0
607
         ind_zmat(izm,5)=0
608
         val_zmat(izm,1)=norm1
609
         val_zmat(izm,2)=0.0
610
         val_zmat(izm,3)=0.0
611
         DejaFait(Izm2)=.TRUE.
612

    
613
         Izm=3
614

    
615
         CALL vecteur(izm1,izm3,x,y,z,vx2,vy2,vz2,norm2)
616
         val=angle(vx1,vy1,vz1,norm1, &
617
              vx2,vy2,vz2,norm2)
618

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

    
621
         ind_zmat(izm,1)=izm3
622
         ind_zmat(izm,2)=izm1
623
         ind_zmat(izm,3)=izm2
624
         ind_zmat(izm,4)=0
625
         ind_zmat(izm,5)=0
626
         val_zmat(izm,1)=norm2
627
         val_zmat(izm,2)=val
628
         val_zmat(izm,3)=0.0
629
         DejaFait(Izm3)=.TRUE.
630

    
631
      END IF
632

    
633
      IF (Lies_CF(Izm1,0).eq.1) THEN
634
         if (debug) WRITE(*,*) 'Cas merdouille,',IdAt0, &
635
              ' li? a 1 atome'
636

    
637
         IZm2=Lies_CF(IdAt0,1)
638
         Izm=2
639
         CALL vecteur(IZm1,IZm2,x,y,z,vx1,vy1,vz1,norm1)
640
!     write(*,11) n1,n2,norm1
641

    
642
         ind_zmat(izm,1)=izm2
643
         ind_zmat(izm,2)=izm1
644
         ind_zmat(izm,3)=0
645
         ind_zmat(izm,4)=0
646
         ind_zmat(izm,5)=0
647
         val_zmat(izm,1)=norm1
648
         val_zmat(izm,2)=0.0
649
         val_zmat(izm,3)=0.0
650
         DejaFait(Izm2)=.TRUE.
651

    
652
!     Pour les autres atomes, on prend les atomes
653
!     CF lie a l'tome CF lie a At0... en esperant
654
!     qu'il en a !!!
655
!     => il faudra voir le cas ou il n'en a pas !!!
656
!     en attendant on rajoute le test...
657
         IF (Lies_CF(Izm2,0).Eq.0) THEN
658
            WRITE(*,*) "Je n'arrive pas a trouver les premiers"
659
            WRITE(*,*) "atomes pour construire la molecule"
660
            WRITE(*,*) "Atome central li? au plus a 1 atome",Izm1,izm2
661
            WRITE(*,*) "Et atome 2 li? a rien... moi perdu"
662
            STOP
663
         END IF
664

    
665
!     On ajoute les atomes lies a cet atome dans la liste a faire
666
         Do iaf=1,Lies_CF(Izm2,0)
667
            CAfaire(IndFin)=Lies_CF(Izm2,Iaf)
668
            IndFin=IndFin+1
669
         END DO
670
         CAfaire(IndFin)=0
671
         Izm3= Lies_CF(Izm2,1)
672
         Izm=3
673
         CALL vecteur(izm2,izm1,x,y,z,vx1,vy1,vz1,norm1)
674

    
675
         CALL vecteur(izm2,izm3,x,y,z,vx2,vy2,vz2,norm2)
676
         val=angle(vx1,vy1,vz1,norm1, &
677
              vx2,vy2,vz2,norm2)
678

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

    
681
         ind_zmat(izm,1)=izm3
682
         ind_zmat(izm,2)=izm2
683
         ind_zmat(izm,3)=izm1
684
         ind_zmat(izm,4)=0
685
         ind_zmat(izm,5)=0
686
         val_zmat(izm,1)=norm1
687
         val_zmat(izm,2)=val
688
         val_zmat(izm,3)=0.0
689
         DejaFait(Izm3)=.TRUE.
690

    
691
!     je pense que ce qui suit est totalement inutile
692
!     C     DO IdAtl=3,Lies_CF(IdAt0,0)
693
!     C     Izm4= Lies_CF(IdAt0,IdAtl)
694
!     C     Izm=Izm+1
695
!     C     Call AddLigne(Izm,Izm4,Izm1,Izm2,Izm3,ind_zmat,val_zmat
696
!     C     $,x,y,z)
697
!     C     DejaFait(Izm4)=.TRUE.
698
!     C     END DO
699
      END IF
700

    
701
!     on a cree la premiere couche autour du premier centre
702
!     reste a finir les autres couches
703
      IndAFaire=1
704
      if (debug) WRITE(IOOUT,*) "CaFaire:",CaFaire
705
      Do WHILE (Cafaire(IndAFaire).ne.0)
706
         IAt0=Cafaire(IndAfaire)
707
         if (debug) WRITE(*,*) "IndAFaire,IAt0,Lies_CF(IAt0,0)", &
708
              IndAFaire,IAt0,Lies_CF(IAt0,0)
709
         if (Lies_CF(IAt0,0).ge.1) THEN
710
!     IAt0 est l'atome pour lequel on construit la couche suivante
711
!     le premier atome lie est particulier car il definit l'orientation
712
!     de ce fragment
713
            IAti=Lies_CF(IAt0,1)
714
!     WRITE(IOOUT,*) 'IAti:',IAti
715
            IAtd=IAt0
716
!     WRITE(IOOUT,*) 'IAtd:',IAtd
717
            IAtv=Lies_CP(IAt0,1)
718
!     WRITE(IOOUT,*) 'IAtv:',IAtv
719
            IIAtdh=1
720
            IAtdh=Lies_CF(IAtv,IIatdh)
721
            DO WHILE (Iatdh.eq.IAt0)
722
               IIatdh=IIatdh+1
723
               IAtdh=Lies_CF(IAtv,IIatdh)
724
            END DO
725
            IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
726
!     WRITE(IOOUT,*) 'IAtdh:',IAtdh
727
            IF (.NOT.DejaFait(IAti)) THEN
728
               Izm=Izm+1
729
               if (debug) WRITE(*,*) 'izm,IAti,IAtd,IAtv,IAtdh', &
730
                    izm,IAti,IAtd,IAtv,IAtdh
731
               Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh,ind_zmat,val_zmat &
732
                    ,x,y,z)
733
               DejaFait(IAti)=.TRUE.
734
            END IF
735

    
736
            CAfaire(IndFin)=Lies_CF(IAt0,1)
737
            IndFin=IndFin+1
738
!     Le traitement des autres est plus facile
739
            IAtdh=Lies_CF(IAt0,1)
740
            DO IAta=2,Lies_CF(IAt0,0)
741
               IAtI=Lies_CF(IAt0,IAta)
742
                     if (debug) WRITE(*,*) " Problem is here; IndFin:",I
743
               CAfaire(IndFin)=IAtI
744
               IndFin=IndFin+1
745

    
746
               IF (.NOT.DejaFait(IAti)) THEN
747
                  Izm=Izm+1
748
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh,ind_zmat &
749
                       ,val_zmat,x,y,z)
750
                  DejaFait(IAti)=.TRUE.
751
               END IF
752
            END DO
753

    
754
            CAfaire(IndFin)=0
755
         END IF
756
         IndAFaire=IndAFaire+1
757
      END DO
758

    
759

    
760
!     On a fini de creer la molecule autour du premier atome 'centre'
761
!     Pour les autres c'est presque id... sauf que les axes sont deja fixes
762
!     On pourrait imaginer de mettre des atomes fictifs pour redefinir des axes
763
!     locaux... dans une autre version
764
!     V2.0
765
      NbAt0=NbAt0-1
766
      DO I=1,NbAt0
767
!     Boucle pour trouver celui qui a le plus d'atomes CF
768

    
769
         IdAt0=0
770
         VCF=0
771

    
772
         DO ICAt=1,na
773
            if ((Lies_CP(ICaT,0).eq.0).AND.(Lies_CF(ICAt,0).gt.VCF)) &
774
                 THEN
775
               IdAt0=ICat
776
               VCF=Lies_CF(ICAt,0)
777
            END IF
778
         END DO
779
         Lies_CP(IdAt0,0)=-1
780

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

    
787
         if (debug) THEN
788
            WRITE(IOOUT,*) 'LIAISONS tralala',nbcycle
789
            DO IAt=1,na
790
               WRITE(IOOUT,1003) iat,(LIAISONS(Iat,Jat),Jat=0,NMaxL)
791
            END DO
792
         END IF
793

    
794
!     Boucle pour voir quel est l'atome du fragment precedent qui est le plus
795
!     proche de celui ci
796
!     Rq: si ce sont des cycles lies, il pourrait etre malin de rechercher
797
!     a quoi il etait lie.
798
         norm2=1.e6
799
         Idx=0
800
         DO ICAt=1,Izm
801
            CALL vecteur(Idat0,ind_zmat(icat,1),x,y,z,vx1,vy1,vz1,norm1)
802
            if (norm2.ge.norm1) THEN
803
               norm2=norm1
804
               Idx=Icat
805
            END IF
806
         END DO
807
         IF (debug) WRITE(*,*) 'Link between fragments ',IdAt0, &
808
              ind_zmat(Idx,1), ' -- Idx=',Idx
809

    
810
!     on a l'indice... on va rajouter cet atome a la liste !
811
         izm=izm+1
812
         n1=ind_zmat(Idx,1)
813
!     Petit probleme si Idx<=2
814
         IF (Idx.EQ.1) THEN
815
!     Pb non negligeable ...
816
            IF (izm.ge.2) THEN
817
               IAtv=ind_zmat(2,1)
818
               IAtdh=ind_zmat(3,1)
819
            ELSE
820
               WRITE(*,*) '2 fragments, le 1er a moins de 2 atomes...'
821
               WRITE(*,*) "J'ai merde... "
822
               STOP
823
            END IF
824
         ELSEIF (Idx.EQ.2) THEN
825
            IAtv=1
826
            IF (izm.ge.2) THEN
827
               IAtdh=ind_zmat(3,1)
828
            ELSE
829
               WRITE(*,*) '2 fragments, le 1er a moins de 2 atomes...'
830
               WRITE(*,*) "J'ai merde... "
831
               STOP
832
            END IF
833
         ELSE
834
            IAtv=ind_zmat(Idx,2)
835
            IAtdh=ind_zmat(Idx,3)
836
         END IF
837

    
838
!     WRITE(*,*) "Outchilou IAtv, IAtdh", IAtv, IAtdh,n1
839
         CALL vecteur(n1,IdAt0,x,y,z,vx1,vy1,vz1,norm1)
840
         CALL vecteur(n1,IAtv,x,y,z,vx2,vy2,vz2,norm2)
841
         val=angle(vx1,vy1,vz1,norm1, &
842
              vx2,vy2,vz2,norm2)
843
         if (debug) WRITE(*,*) 'Angle val:1',IdAt0,n1,IAtv,val
844
         if ((abs(val).LE.10.).OR.(abs(pi-val).le.10.)) THEN
845
            IAtv=IAtdh
846
!     Ceci pose pb si Idx<=3... a traiter plus tard
847
            IF (Idx.ge.3) THEN
848
               IAtdh=ind_zmat(Idx,4)
849
            ELSE
850
               if (izm.ge.4) THEN
851
                  IAtdh=4
852
               ELSE
853
                  WRITE(*,*) 'Pb Angle Val+idx=3+izm<=4'
854
                  STOP
855
               END IF
856
            END IF
857
         END IF
858
         Call AddLigne(Izm,IdAt0,ind_zmat(Idx,1),IAtv,IAtdh, &
859
              ind_zmat,val_zmat,x,y,z)
860

    
861

    
862
         IndFin=1
863
         IAtd=IdAt0
864
         n1=IAtdh
865
         IAtdh=IAtv
866
         IAtv=ind_zmat(Idx,1)
867
!     On ajoute les atomes lies a cet atome dans la liste a faire
868
         Do iaf=1,Lies_CF(IdAt0,0)
869
            IatI=Lies_CF(IdAt0,Iaf)
870
!     We check that this atom is not the one that is the closest to
871
!     the center atom
872
            if (IAtv.Ne.IAtI) THEN
873
               IF (debug) WRITE(*,*) 'Adding atom',IAtI, &
874
                    'to CAFaire in pos',iaf
875
               CAfaire(IndFin)=IAtI
876
               IndFin=IndFin+1
877
!     On les ajoute aussi dans la zmat
878
               If (.NOT.DejaFait(IatI)) THEN
879
                  izm=izm+1
880
!     WRITE(*,*) "Outchili IAtv, IAtdh,IAtI", IAtv, IAtdh,IAtI
881
                  CALL vecteur(IAtI,IAtD,x,y,z,vx1,vy1,vz1,norm1)
882
                  CALL vecteur(IAtI,IAtv,x,y,z,vx2,vy2,vz2,norm2)
883
                  val=angle(vx1,vy1,vz1,norm1, vx2,vy2,vz2,norm2)
884
                  if (debug) WRITE(*,*) 'Angle val:2',IAtI,IAtD,IAtv,val
885
                  if ((abs(val).LE.10.).OR.(abs(180.-abs(val)).le.10.)) &
886
                       THEN
887
                     Call AddLigne(Izm,IAtI,IAtv,IAtdh,n1, &
888
                          ind_zmat,val_zmat,x,y,z)
889
                  ELSE
890

    
891
                     Call AddLigne(Izm,IAtI,IAtD,IAtv,IAtdh, &
892
                          ind_zmat,val_zmat,x,y,z)
893
                  END IF
894
               END IF
895
            END IF
896
         END DO
897
         CAfaire(IndFin)=0
898

    
899
!     On traite le reste de ce fragment !!
900
         IndAFaire=1
901
         WRITE(IOOUT,*) CaFaire
902
         Do WHILE (Cafaire(IndAFaire).ne.0)
903
            IAt0=Cafaire(IndAfaire)
904
            if (Lies_CF(IAt0,0).ge.1) THEN
905
!     IAt0 est l'atome pour lequel on construit la couche suivante
906
!     le premier atome lie est particulier car il definit l'orientation
907
!     de ce fragment
908
               Itmp=1
909
               IAti=Lies_CF(IAt0,Itmp)
910
               DO WHILE (DejaFait(IAti).AND.(Itmp.Le.Lies_CF(IAt0,0)))
911
                  Itmp=Itmp+1
912
                  IAti=Lies_CF(IAt0,ITmp)
913
               END DO
914
               If (.NOT.DejaFait(Iati)) THEN
915
                  IF (debug) WRITE(IOOUT,*) 'IAti:',IAti
916
                  IAtd=IAt0
917
                  IF (debug)     WRITE(IOOUT,*) 'IAtd:',IAtd
918
                  IAtv=Lies_CP(IAt0,1)
919
                  IF (debug)     WRITE(IOOUT,*) 'IAtv:',IAtv
920
                  IIAtdh=1
921
                  IAtdh=Lies_CF(IAtv,IIatdh)
922
                  DO WHILE (Iatdh.eq.IAt0)
923
                     IIatdh=IIatdh+1
924
                     IAtdh=Lies_CF(IAtv,IIatdh)
925
                  END DO
926
                  IF (IAtdh.eq.0) IAtdh=Lies_CP(IAtv,1)
927
                  IF (debug) WRITE(IOOUT,*) 'IAtdh:',IAtdh
928
                  Izm=Izm+1
929
                  Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
930
                       ind_zmat,val_zmat,x,y,z)
931
!     Le traitement des autres est plus facile
932
                  IAtdh=Lies_CF(IAt0,ITmp)
933
                  DO IAta=ITmp+1,Lies_CF(IAt0,0)
934
                     IAtI=Lies_CF(IAt0,IAta)
935
                     If (.NOT.DejaFait(IAtI)) THEN
936
                        Izm=Izm+1
937
                        Call AddLigne(Izm,IAti,IAtd,IAtv,IAtdh, &
938
                             ind_zmat,val_zmat,x,y,z)
939
                     END IF
940
                  END DO
941
               END IF
942

    
943
!     On ajoute les atomes lies a cet atome dans la liste a faire
944
               Do iaf=1,Lies_CF(IAt0,0)
945
                  CAfaire(IndFin)=Lies_CF(IAt0,Iaf)
946
                  IndFin=IndFin+1
947
               END DO
948
               CAfaire(IndFin)=0
949
            END IF
950
            IndAFaire=IndAFaire+1
951
         END DO
952
!12345    CONTINUE
953
      END DO
954

    
955
      if (debug) THEN
956
         WRITE(*,*) 'ind_zmat'
957
         DO I=1,na
958
            WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
959
         END DO
960
      END IF
961

    
962

    
963

    
964
      if (debug) WRITE(*,*) "====================== Exiting Calc_zmat ==================="
965

    
966
    END SUBROUTINE Calc_Zmat
967