Statistiques
| Révision :

root / src / Calc_zmat_constr_frag.f90 @ 8

Historique | Voir | Annoter | Télécharger (66,62 ko)

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

    
4
  !     
5
  !     This second version enables the user to freeze some atoms
6
  !     the frozen atoms indices are listed in the frozen array.
7
  !     
8
  !     The idea is surely too stupid to work all the time... but here it is
9
  !     we first construct the zmat using only the frozen atoms, and then
10
  !     we add the other atoms on top of the first ones...
11
  !     
12
  Use Path_module, only : max_Z, NMaxL, Nom
13
  Use Io_module
14

    
15
  IMPLICIT NONE
16

    
17

    
18
  CHARACTER(5) :: AtName
19
  integer(KINT) :: na,atome(na),ind_zmat(Na,5)
20
  INTEGER(KINT) :: idx_zmat(NA)
21
  real(KREAL) ::  x(Na),y(Na),z(Na),fact
22
  real(KREAL) ::  val_zmat(Na,3)
23
  real(KREAL) :: r_cov(0:Max_Z)
24
  !     Frozen contains the indices of frozen atoms
25
  INTEGER(KINT) :: Frozen(*),NFroz
26
  INTEGER(KINT), ALLOCATABLE :: FrozDist(:) ! na
27
  INTEGER(KINT) :: NbFrag,IdxFrag
28
  INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !na
29
  INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !na,na
30
  INTEGER(KINT), ALLOCATABLE :: FrozFrag(:,:) !na,3
31
  !  INTEGER(KINT), ALLOCATABLE :: IdxFragAt(:) !na
32
  INTEGER(KINT), ALLOCATABLE :: FrozBlock(:,:) !(na,0:na)
33
  REAL(KREAL), ALLOCATABLE :: DistFroz(:) !na
34

    
35
  INTEGER(KINT) :: IdxCaFaire, IAfaire
36
  INTEGER(KINT), ALLOCATABLE ::  LIAISONS(:,:) ! (Na,0:NMaxL)
37
  INTEGER(KINT), ALLOCATABLE ::  LiaisonsBis(:,:) ! (Na,0:NMaxL)
38
  INTEGER(KINT), ALLOCATABLE ::  LiaisonsIni(:,:) ! (Na,0:NMaxL)
39
  INTEGER(KINT), ALLOCATABLE :: CAFaire(:) ! (Na)
40

    
41

    
42
  real(KREAL) ::  vx1,vy1,vz1,norm1
43
  real(KREAL) ::  vx2,vy2,vz2,norm2
44
  real(KREAL) ::  vx3,vy3,vz3,norm3
45
  real(KREAL) ::  vx4,vy4,vz4,norm4
46
  real(KREAL) ::  vx5,vy5,vz5,norm5
47
  real(KREAL) ::  val,val_d, d12, d13,d23,d
48
  Logical Debug
49
  LOGICAL, ALLOCATABLE :: DejaFait(:), FCaf(:) !(na)
50
  Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen
51
  LOGICAL F1213, F1223,F1323
52

    
53

    
54
  INTEGER(KINT) :: I,J,n0,n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax
55
  INTEGER(KINT) :: IMax, I3,I1,Ip, IFragFroz
56
  INTEGER(KINT) :: I0, Izm, JAt,Izm0
57

    
58
  REAL(KREAL) :: sDihe, Pi
59

    
60
  INTERFACE
61
     function valid(string) result (isValid)
62
       CHARACTER(*), intent(in) :: string
63
       logical                  :: isValid
64
     END function VALID
65

    
66
     FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
67
       use Path_module, only :  Pi,KINT, KREAL
68
       real(KREAL) ::  v1x,v1y,v1z,norm1
69
       real(KREAL) ::  v2x,v2y,v2z,norm2
70
       real(KREAL) ::  angle
71
     END FUNCTION ANGLE
72

    
73
     FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
74
       use Path_module, only :  Pi,KINT, KREAL
75
       real(KREAL) ::  v1x,v1y,v1z,norm1
76
       real(KREAL) ::  v2x,v2y,v2z,norm2
77
       real(KREAL) ::  SinAngle
78
     END FUNCTION SINANGLE
79

    
80

    
81
     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
82
       use Path_module, only :  Pi,KINT, KREAL
83
       real(KREAL) ::  v1x,v1y,v1z,norm1
84
       real(KREAL) ::  v2x,v2y,v2z,norm2
85
       real(KREAL) ::  v3x,v3y,v3z,norm3
86
       real(KREAL) ::  angle_d,ca,sa
87
     END FUNCTION ANGLE_D
88

    
89
  END INTERFACE
90

    
91

    
92

    
93
  Pi=dacos(-1.d0)
94
  debug=valid("calc_zmat_constr").OR.valid("calc_zmat_contr_frag")
95

    
96
  if (debug) WRITE(*,*) "=============================== Entering Calc_zmat_contr_frag ========================"
97
  if (na.le.2) THEN
98
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
99
     RETURN
100
  END IF
101

    
102
  ALLOCATE(FrozDist(Na),Fragment(na), NbAtFrag(na),FragAt(na,na))
103
  !  ALLOCATE(FrozFrag(na,3), IdxFragAt(na), FrozBlock(na,0:na))
104
  ALLOCATE(FrozFrag(na,3), FrozBlock(na,0:na))
105
  ALLOCATE(DistFroz(na),Liaisons(na,0:NMaxL))
106
  ALLOCATE(LiaisonsBis(na,0:NMaxL),LiaisonsIni(na,0:NMaxL))
107
  ALLOCATE(CaFaire(na),DejaFait(Na),FCaf(na),FrozAt(na))
108

    
109
  DO i=1,na
110
     DO J=1,5
111
        Ind_Zmat(i,J)=0
112
     END DO
113
  END DO
114

    
115
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
116
!
117
!     Easy case : 3 or less atoms
118
!
119
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
120
  if (Na.eq.3) THEN
121
     d12=sqrt((x(1)-x(2))**2+(y(1)-y(2))**2+(z(1)-z(2))**2)
122
     d13=sqrt((x(1)-x(3))**2+(y(1)-y(3))**2+(z(1)-z(3))**2)
123
     d23=sqrt((x(3)-x(2))**2+(y(3)-y(2))**2+(z(3)-z(2))**2)
124
     F1213=(d12<=d13)
125
     F1323=(d13<=d23)
126
     F1223=(d12<=d23)
127
     if (debug) THEN
128
        WRITE(*,*) "DEBUG Calc_Zmat 3 atoms"
129
        WRITE(*,*) "d12,d13,d23:",d12,d13,d23
130
        WRITE(*,*) "F1213,F1323,F1223:",F1213,F1323,F1223
131
     END IF
132
     if (F1213) THEN
133
        if (F1323) THEN
134
           ind_zmat(1,1)=3
135
           ind_zmat(2,1)=1
136
           ind_zmat(2,2)=3
137
           ind_zmat(3,1)=2
138
           ind_zmat(3,2)=1
139
           ind_zmat(3,3)=3
140
        ELSE
141
           ind_zmat(1,1)=1
142
           ind_zmat(2,1)=2
143
           ind_zmat(2,2)=1
144
           ind_zmat(3,1)=3
145
           ind_zmat(3,2)=2
146
           ind_zmat(3,3)=1
147
        END IF
148
     ELSE
149
        IF (F1223) THEN
150
           ind_zmat(1,1)=2
151
           ind_zmat(2,1)=1
152
           ind_zmat(2,2)=2
153
           ind_zmat(3,1)=3
154
           ind_zmat(3,2)=1
155
           ind_zmat(3,3)=2
156
        ELSE
157
           ind_zmat(1,1)=2
158
           ind_zmat(2,1)=3
159
           ind_zmat(2,2)=2
160
           ind_zmat(3,1)=1
161
           ind_zmat(3,2)=3
162
           ind_zmat(3,3)=2
163
        END IF
164
     END IF
165
     IF (debug) THEN
166
        DO i=1,na
167
           WRITE(*,'(1X,5(1X,I5))') (ind_zmat(i,j),j=1,5)
168
        END DO
169
     END IF
170

    
171
     !     We have ind_zmat, we fill val_zmat
172
     val_zmat(1,1)=0.
173
     val_zmat(1,2)=0.
174
     val_zmat(1,3)=0.
175
     n2=ind_zmat(2,1)
176
     n1=ind_zmat(2,2)
177
     d=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2+(z(n1)-z(n2))**2)
178
     val_zmat(2,1)=d
179
     val_zmat(2,2)=0.
180
     val_zmat(2,3)=0.
181
     n1=ind_zmat(3,1)
182
     n2=ind_zmat(3,2)
183
     n3=ind_zmat(3,3)
184
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
185
     if (debug) write(*,*) n1,n2,norm1
186
     val_zmat(3,1)=norm1
187

    
188
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
189
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
190
     if (debug) write(*,*) n2,n3,norm2,val
191
     val_zmat(3,2)=val
192
     val_zmat(3,3)=0.
193

    
194
     RETURN
195
  END IF
196
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
197
!
198
!  End of   Easy case : 3 or less atoms
199
!
200
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
201
!
202
! General Case 
203
!
204
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
205
!
206
! 1 - Frozen Atoms
207

    
208
! Initialization
209
  DejaFait=.False.
210
  FrozAt=.False.
211
  Liaisons=0
212
  LiaisonsBis=0
213
  ind_zmat=0
214
  val_zmat=0.d0
215

    
216
  NFroz=0
217
  I=1
218
  DO WHILE (Frozen(I).NE.0)
219
     If (Frozen(I).LT.0) THEN
220
        DO J=Frozen(I-1),abs(Frozen(I))
221
           IF (.NOT.FrozAt(J)) THEN
222
              NFroz=NFroz+1
223
              FrozAt(J)=.TRUE.
224
           END IF
225
        END DO
226
     ELSEIF (.NOT.FrozAt(Frozen(I))) THEN
227
        FrozAt(Frozen(I))=.TRUE.
228
        NFroz=NFroz+1
229
     END IF
230
  I=I+1
231
  END DO
232

    
233
  if (debug) WRITE(*,*) 'DGB zmtc NFroz=', NFroz
234
  if (debug) WRITE(*,*) 'DGB zmtc FrozAt=',(FrozAt(I),I=1,na)
235

    
236
  if (debug) THEN
237
     WRITE(*,*) "Liaisons initialized"
238
     WRITE(*,*) 'Covalent radii used'
239
     DO iat=1,na
240
        i=atome(iat)
241
        WRITE(*,*) Nom(I),Iat,r_cov(i),r_cov(i)*fact
242
     END DO
243
  END IF
244

    
245
1003 FORMAT(1X,I4,' - ',25(I5))
246
  !     DO IL=1,na
247
  !     WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
248
  !     WRITE(IOOUT,1003) Il,(LIAISONSBis(IL,JL),JL=0,NMaxL)
249
  !     END DO
250
!   First step : connectivity  are calculated
251

    
252
  Call CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)
253

    
254
  if (debug) THEN
255
     WRITE(*,*) "Connections calculated"
256
     DO IL=1,na
257
        WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
258
     END DO
259
  END IF
260

    
261
  FCaf=.TRUE.
262
  Call Decomp_frag(na,liaisons,FCaf,nbfrag,Fragment,NbAtFrag,FragAt)
263

    
264
  IF (debug) THEN
265
     WRITE(*,*) 'I found ',NbFrag, 'fragments'
266
     WRITE(*,'(20(1X,I4))') (NbAtFrag(I),I=1,NbFrag)
267
     DO I=1,NbFrag
268
        WRITE(*,*) NbAtFrag(I)
269
        WRITE(*,*) 'Fragment ', i
270
        DO J=1,Na
271
           IF (Fragment(J).EQ.I)   THEN                
272
            if (FrozAt(J))    THEN   
273
                   WRITE(*,'(1X,I3,"f",1X,A5,3(1X,F10.6))') J,Nom(Atome(J)) &
274
                        ,X(J),Y(J),Z(J)
275
               ELSE
276
                    WRITE(*,'(1X,I3,2X,A5,3(1X,F10.6))') J,Nom(Atome(J)) &
277
                      ,X(J),Y(J),Z(J)
278
                END IF
279
       END IF
280
        END DO
281
        WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
282
     END DO
283
  END IF
284

    
285

    
286
! First, we decompose the connectivity into Frozen atoms and non frozen atoms
287
  DO I=1,na
288
     LiaisonsBis(I,0)=0
289
     DO J=1,Liaisons(I,0)
290
        IF (FrozAt(Liaisons(I,J))) THEN
291
           LiaisonsBis(I,0)=LiaisonsBis(I,0)+1
292
           LiaisonsBis(I,LiaisonsBis(I,0))=Liaisons(I,J)
293
        END IF
294
     END DO
295
     IMax=LiaisonsBis(I,0)+1
296
     LiaisonsBis(I,Imax)=0
297
     DO J=1,Liaisons(I,0)
298
        IF (.NOT.FrozAt(Liaisons(I,J))) THEN
299
           LiaisonsBis(I,IMax)=LiaisonsBis(I,Imax)+1
300
           LiaisonsBis(I,LiaisonsBis(I,Imax)+Imax)=Liaisons(I,J)
301
        END IF
302
     END DO
303
  END DO
304

    
305
  if (debug) THEN
306
     WRITE(*,*) "Liaisons and LiaisonsBis"
307
     DO I=1,Na
308
        WRITE(*,'(1X,I3," -",I3,":",15(1X,I3))') I, &
309
             (Liaisons(I,J),J=0,NMaxL)
310
        WRITE(*,'(1X,I3," +",I3,":",15(1X,I3))') I, &
311
             (LiaisonsBis(I,J),J=0,NMaxL)
312
     END DO
313
  END IF
314

    
315
! Now, for each frozen atom, we count the length of the frozen block
316
! FrozBlock(I,0) contains the number of frozen atoms forming a frozen
317
!                fragment containing atom I
318
! FrozBlock(I,:) is the list of the frozen atoms of this fragment.
319
  Do I=1,na
320
     FrozBlock(I,0)=-1
321
  END DO
322
  DO I=1,na
323
     IF (FrozAt(I).AND.(FrozBlock(I,0).EQ.-1)) THEN
324
        if (debug) WRITE(*,*) 'Treating atom ',I
325
        FrozBlock(I,0)=1
326
        FrozBlock(I,1)=I
327
        DO J=1,na
328
           DejaFait(J)=.FALSE.
329
        END DO
330
        DejaFait(I)=.TRUE.
331
        DO J=1,LiaisonsBis(I,0)
332
           CaFaire(J)=LiaisonsBis(I,J)
333
        END DO
334
        IdxCaFaire=LiaisonsBis(I,0)+1
335
        CaFaire(IdxCaFaire)=0
336
        IAfaire=1
337
        FCaf=DejaFait
338
        DO WHILE (CaFaire(IAfaire).NE.0)
339
           Iat=CaFaire(IAFAire)
340
           if (debug) WRITE(*,*) 'IaFaire; Iat:', IaFaire, Iat, DejaFait(Iat)
341
           IF (.NOT.DejaFait(Iat)) THEN
342
              FrozBlock(I,0)=FrozBlock(I,0)+1
343
              FrozBlock(I,FrozBlock(I,0))=Iat
344
              DO J=1,LiaisonsBis(Iat,0)
345
                 IF ((.NOT.DejaFait(LiaisonsBis(Iat,J))).AND.(.NOT.FCaf(LiaisonsBis(Iat,J)))) THEN
346
                    CaFaire(IdxCaFaire)=LiaisonsBis(Iat,J)
347
                    IdxCaFaire=IdxCaFaire+1
348
                    CaFaire(IdxCaFaire)=0
349
                    FCaf(LiaisonsBis(Iat,J))=.TRUE.
350
                 END IF
351
              END DO
352
              !     WRITE(*,*) 'liaisonbis(Iat,0),FrozBlick(I,0)',&
353
              !                      LiaisonsBis(Iat,0), FrozBlock(I,0)
354
           END IF
355
           DejaFait(Iat)=.TRUE.
356
           IaFaire=IaFaire+1
357
        END DO
358
        if (debug) WRITE(*,*) 'I,FrozBlock(0),IaFaire',I,FrozBlock(I,0),IaFaire
359
        if (debug) WRITE(*,*) 'FrozBlock:',FrozBlock(I,1:FrozBlock(I,0))
360
        !        FrozBlock(I,1)=I
361
        !        DO J=2,IaFaire
362
        !           FrozBlock(I,J)=CaFaire(J-1)
363
        !        END DO
364
        !     We copy this block into all frozen atoms that belongs to it
365
        DO J=2,Frozblock(I,0)
366
           Iat=FrozBlock(I,J)
367
           DO K=0,FrozBlock(I,0)
368
              FrozBlock(Iat,K)=FrozBlock(I,K)
369
           END DO
370
        END DO
371
     ELSE
372
        IF (.NOT.FrozAt(I))      FrozBlock(I,0)=0
373
     END IF
374
  END DO
375

    
376
  if (debug) THEN
377
     WRITE(*,*) "FrozBlock"
378
     DO I=1,Na
379
        If (FrozAt(I))  WRITE(*,'(1X,I3," -",I3,":",15(1X,I3))') I, &
380
             (FrozBlock(I,J),J=0,FrozBlock(I,0))
381
     END DO
382
  END IF
383

    
384
  DO I=1,NbFrag
385
     FrozFrag(I,1)=0
386
     FrozFrag(I,2)=0
387
     DO J=1,NbAtFrag(I)
388
        IF (FrozAt(FragAt(I,J))) THEN
389
           FrozFrag(I,1)=FrozFrag(I,1)+1
390
           IF (FrozBlock(FragAt(I,J),0).GE.FrozFrag(I,2))  &
391
                FrozFrag(I,2)=FrozBlock(FragAt(I,J),0)
392
           if (FrozFrag(I,3).LE.LiaisonsBis(FragAt(I,J),0))& 
393
                FrozFrag(I,3)=LiaisonsBis(FragAt(I,J),0)
394
        END IF
395
     END DO
396
     IF (debug) WRITE(*,*) 'Frag :',I,FrozFrag(I,1), &
397
          ' frozen atoms,max linked:'  &
398
          ,FrozFrag(I,2),' with at most',FrozFrag(I,3),' frozen'
399
  END DO
400

    
401
  !     We will now build the molecule fragment by fragment
402
  !     First the frozen atoms, then the rest of the fragment
403
  !     We choose the starting fragment with two criteria:
404
  !     1- Number of linked frozen atoms:
405
  !       * >=3 is good as it fully defines the coordinate space
406
  !       * 2 is ok as we can either use a 3rd atom from the same fragment
407
  !       or add a X atom somewhere but this complicates quite a lot the way
408
  !       to treat the conversion from cartesian to zmat latter
409
  !       * 1 is bad...
410
  !     2- Max number of linked frozen atoms
411
  !     this allows us to deal more easily with cases 1- when number of 
412
  !     directly linked frozen atoms is less than 3
413

    
414
  IFrag=0
415
  I0=0
416
  I1=0
417
  DO I=1,NbFrag
418
     SELECT CASE(FrozFrag(I,3)-I0)
419
     CASE (1:)
420
        IFrag=I
421
        I0=FrozFrag(I,3)
422
        I1=FrozFrag(I,2)
423
     CASE (0)
424
        if (FrozFrag(I,2).GT.I1) THEN
425
           IFrag=I
426
           I0=FrozFrag(I,3)
427
           I1=FrozFrag(I,2)
428
        END IF
429
     END SELECT
430
  END DO
431

    
432
  if (debug) WRITE(*,'(1X,A,I5,A,I5,A,I5,A)') 'Starting with fragment:',IFrag,' with ',I0,' frozen linked and overall',I1,' linked'
433

    
434
  !     We will build the first fragment in a special way, as it will
435
  !     set the coordinates system
436
  !     We look for the frozen atom that is linked to the maximum frozen atom,
437
  !     and belongs to the longer fragment
438
  I0=0
439
  I1=0
440
  DO I=1,NbAtFrag(IFrag)
441
     IF (FrozAt(FragAt(IFrag,I))) THEN
442
        SELECT CASE(LiaisonsBis(FragAt(IFrag,I),0)-I0)
443
        CASE (1:) 
444
           IAt=FragAt(IFrag,I)
445
           I0=LiaisonsBis(IAt,0)
446
           I1=FrozBlock(IAt,0)
447
        CASE (0)
448
           IF (FrozBlock(FragAt(IFrag,I),0).GT.I1) THEN
449
              IAt=FragAt(IFrag,I)
450
              I0=LiaisonsBis(Iat,0)
451
              I1=FrozBlock(Iat,0)
452
           END IF
453
        END SELECT
454
     END IF
455
     if (debug) WRITE(*,*) 'DBG: I,FragAt(IFrag,I),IAt,I0,I1',I,FragAt(IFrag,I),IAt,I0,I1
456
  END DO
457

    
458
  if (debug) WRITE(*,*) 'Fragment 1, starting with atom:',IAt,I0,I1
459

    
460
  DO I=1,na
461
     DejaFait(I)=.FALSE.
462
     FCaf(I)=.FALSE.
463
  END DO
464

    
465
  izm=0
466
  SELECT CASE (I0)
467
  CASE(3:)
468
     if (debug) WRITE(*,*) 'DBG select case I0 3'
469
     n0=Iat
470
     !     We search for the fourth atom while making sure that the dihedral angle
471
     !     is not 0 or Pi
472
     ITmp=2
473
     sDihe=0.
474
     n2=IAt
475
     n3=LiaisonsBis(Iat,1)
476
     ! We search for the third atom while making sure that it is not aligned with first two !
477
     DO While ((ITmp.LE.LiaisonsBis(Iat,0)).AND.(sDihe.LE.0.09d0))
478
        n4=LiaisonsBis(Iat,Itmp)
479
        CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
480
        CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
481
        val_d=angle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
482
        sDiHe=abs(sin(val_d*pi/180.d0))
483
        if (debug) Write(*,*) 'Trying n2,n3,n4,sdihe,val_d',n2,n3,n4,sdihe,val_d
484
        Itmp=Itmp+1
485
     END DO
486
     If (debug) WRITE(*,*) 'Itmp,LiaisonsBis',Itmp,LiaisonsBis(Iat,1:NMaxL)
487
     LiaisonsBis(Iat,Itmp-1)=LiaisonsBis(iat,2)
488
     LiaisonsBis(Iat,2)=n4
489
     If (debug) WRITE(*,*) 'Itmp,LiaisonsBis',Itmp,LiaisonsBis(Iat,1:NMaxL)
490

    
491
     if (sDihe.LE.0.09d0) THEN
492
        WRITE(*,*) "Dans le caca car tous les atoms de ce block sont align?s!"
493
     END IF
494

    
495
     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
496
          vx5,vy5,vz5,norm5)
497

    
498

    
499
     Itmp=2
500
     sDiHe=0.
501

    
502
     DO While ((ITmp.LT.LiaisonsBis(Iat,0)).AND.(sDihe.LE.0.09d0))
503
        ITmp=ITmp+1
504
        n1=LiaisonsBis(Iat,Itmp)
505
        if (debug) WRITe(*,*) 'trying n1,n2,n3,n4',n1,n2,n3,n4
506
        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
507
        ! Is this atom aligned with n2-n3 ?
508
        val_d=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
509
        sDiHe=abs(sin(val_d*pi/180.d0))
510
        if (debug) Write(*,*) 'Angle n3-n2-n1',val_d
511
        if (sDiHe.le.0.09d0) THEN
512
           ! As atoms n2,n3 and n4 are not aligned, we interchange n3 and n4 so that n1,n2 and n3 are not aligned
513
           if (debug) WRITE(*,*) "n3-n2-n1 aligned, we interchange n3 and n4"
514
           n1=n3
515
           n3=n4
516
           n4=n1
517
           n1=LiaisonsBis(Iat,ITmp)
518
           CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)           
519
           CALL vecteur(n2,n4,x,y,z,vx3,vy3,vz3,norm3)           
520
           val_d=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
521
           if (debug) Write(*,*) 'NEW Angle n3-n2-n1',val_d
522
        END IF
523

    
524
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
525
        ! avant de tester l'angle diedre, il faut verifier que ce 4e atome n'est pas
526
        ! aligne avec les 2 premiers. 
527
        ! comme on a deja test? que les 3 premiers ne sont pas alignes,
528
        ! s'il est align? avec les 2 premiers, on peut inverser le role de 2 et 3.
529
        ! On pourrait tout simplier en faisant une bete recherche parmi tous les atomes geles
530
        ! de ce bloc (au moins 4 ?) avec le critere 1) on les range par distance croissante
531
        ! 2) on les scanne tant que l'angle valence n'est pas bon, puis tant que diehedre pas bon
532
        ! on pourrait comme ca faire un tableau avec les atomes ranges d'abord pour le fragment s?lectionn?
533
        ! puis les atomes des autres fragment par distance croissante.
534
        ! les autres fragments ne seraient additonn?s que si l'on ne trouve pas notre bonheur dans le premier bloc
535
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
536

    
537
        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
538
             vx4,vy4,vz4,norm4)
539
        val_d=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
540
             vx2,vy2,vz2,norm2)
541
        sDihe=abs(sin(val_d*pi/180.d0))
542
        if (debug) WRITE(*,*) 'n2,n3,n4,n1, angle_d',n2,n3,n4,n1,val_d
543
     END DO
544

    
545
     DejaFait(n2)=.TRUE.
546
     DejaFait(n3)=.TRUE.
547
     DejaFait(n4)=.TRUE.
548

    
549
     if (sDihe.LE.0.09d0) THEN
550
        !     None of the frozen atoms linked to IAt are good to define the third
551
        !     direction in space...
552
        !     We will look at the other frozen atoms
553
        !     we might improve the search so as to take the atom closest to IAt
554
        if (debug) WRITE(*,*) "All atoms linked to ",Iat," are in a plane. Looking for other frozen atoms"
555
        ITmp=0
556
        DO I=1,NbAtFrag(IFrag)
557
           JAt=FragAt(IFrag,I)
558
           if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
559
              n1=JAt
560
              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
561
              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
562
                   vx4,vy4,vz4,norm4)
563
              val_d=angle_d(vx4,vy4,vz4,norm4, &
564
                   vx5,vy5,vz5,norm5, &
565
                   vx2,vy2,vz2,norm2)
566
              if (abs(sin(val_d)).GE.0.09D0) THEN
567
                 ITmp=ITmp+1
568
                 DistFroz(ITmp)=Norm1
569
                 FrozDist(ITmp)=JAt
570
              END IF
571
           END IF
572
        END DO
573
        IF (ITmp.EQ.0) THEN
574
           !     All dihedral angles between frozen atoms are less than 5?
575
           !     (or more than 175?). We have to look at other fragments :-(
576
           DO I=1,NFroz
577
              Jat=Frozen(I)
578
              if (.NOT.DejaFait(JAt)) THEN
579
                 n1=JAt
580
                 CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
581
                 CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
582
                      vx4,vy4,vz4,norm4)
583
                 val_d=angle_d(vx4,vy4,vz4,norm4, &
584
                      vx5,vy5,vz5,norm5, &
585
                      vx2,vy2,vz2,norm2)
586
                 if (abs(sin(val_d)).GE.0.09D0) THEN
587
                    ITmp=ITmp+1
588
                    DistFroz(ITmp)=Norm1
589
                    FrozDist(ITmp)=JAt
590
                 END IF
591
              END IF
592
           END DO
593
           IF (ITmp.EQ.0) THEN
594
              !     All frozen atoms are in a plane... too bad
595
              WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
596
              WRITE(*,*) 'For now, I do not treat this case'
597
              STOP
598
           END IF
599
        END IF
600
        I1=0
601
        d=1e9
602
        DO I=1,ITmp
603
           IF (DistFroz(I).LE.d) THEN
604
              d=DistFroz(I)
605
              I1=FrozDist(I)
606
           END IF
607
        END DO
608
     ELSE                !(sDihe.LE.0.09d0)
609
        I1=FrozBlock(IAt,ITmp)
610
        if (debug) WRITE(*,*) 'I1,n1:',I1,n1
611
     END IF                 !(sDihe.LE.0.09d0)
612
     !     we now have the atom that is closer to the first one and that
613
     !     forms a non 0/Pi dihedral angle
614

    
615
     ind_zmat(1,1)=n2
616
     ind_zmat(2,1)=n3
617
     ind_zmat(2,2)=n2
618
     ind_zmat(3,1)=n4
619
     ind_zmat(3,2)=n2
620
     ind_zmat(3,3)=n3
621
     DejaFait(n2)=.TRUE.
622
     DejaFait(n3)=.TRUE.
623
     DejaFait(n4)=.TRUE.
624
     CaFaire(1)=n3
625
     CaFaire(2)=n4
626

    
627
     ind_zmat(4,1)=I1
628
     ind_zmat(4,2)=n2
629
     ind_zmat(4,3)=n3
630
     ind_zmat(4,4)=n4
631
     DejaFait(I1)=.TRUE.
632
     CaFaire(3)=I1
633
     CaFaire(4)=0
634
     IdxCaFaire=4
635

    
636
     FCaf(n2)=.TRUE.
637
     FCaf(n3)=.TRUE.
638
     FCaf(I1)=.TRUE.
639
     izm=4
640
     DO I=3,LiaisonsBis(Iat,0)
641
        IF (.NOT.DejaFait(LiaisonsBis(Iat,I))) THEN
642
           izm=izm+1
643
           !           ind_zmat(izm,1)=LiaisonsBis(Iat,I)
644
           !           ind_zmat(izm,2)=n2
645
           !           ind_zmat(izm,3)=n3
646
           !           ind_zmat(izm,4)=n4
647
           Call add2indzmat(na,izm,LiaisonsBis(Iat,I),n2,n3,n4,ind_zmat,x,y,z)
648
           IF (.NOT.FCaf(LiaisonsBis(Iat,I))) THEN
649
              CaFaire(IdxCaFaire)=LiaisonsBis(Iat,I)
650
              IdxCaFaire=IdxCaFaire+1
651
              CaFaire(IdxCaFaire)=0
652
              FCaf(LiaisonsBis(Iat,I))=.TRUE.
653
           END IF
654
           DejaFait(LiaisonsBis(Iat,I))=.TRUE.
655
        END IF
656
     END DO
657

    
658
     if (debug) THEN
659
        WRITE(*,*) "Ind_zmat 0 - SELECT CASE I0 3: -- izm=",izm
660
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
661
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
662
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
663
        DO I=4,izm
664
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2),  &
665
                ind_zmat(I,3), ind_zmat(I,4)
666
        END DO
667
     END IF
668

    
669

    
670
     !     First four atoms (at least) have been put we go on with next parts
671
     !     of this fragment... later
672

    
673

    
674
  CASE(2)
675
     if (debug) WRITE(*,*) 'DBG select case I0 2'
676
     if (debug) WRITE(*,*) 'ATTENTION For now, case with only 3 frozen atoms non treated'
677
     ind_zmat(1,1)=IAt
678
     ind_zmat(2,1)=liaisonsBis(IAt,1)
679
     ind_zmat(2,2)=IAt
680
     ind_zmat(3,1)=LiaisonsBis(IAt,2)
681
     ind_zmat(3,2)=IAt
682
     ind_zmat(3,3)=LiaisonsBis(IAt,1)
683
     DejaFait(IAt)=.TRUE.
684
     DejaFait(LiaisonsBis(Iat,1))=.TRUE.
685
     DejaFait(LiaisonsBis(Iat,2))=.TRUE.
686
     CaFaire(1)=LiaisonsBis(Iat,1)
687
     CaFaire(2)=LiaisonsBis(Iat,2)
688
     FCaf(LiaisonsBis(Iat,1))=.TRUE.
689
     FCaf(LiaisonsBis(Iat,2))=.TRUE.
690

    
691
     !     We search for a fourth atom, first in the FrozBlock atoms
692
     ITmp=2
693
     sDihe=0.
694
     n2=IAt
695
     n3=LiaisonsBis(Iat,1)
696
     n4=LiaisonsBis(Iat,2)
697
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
698
     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
699
     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
700
          vx5,vy5,vz5,norm5)
701

    
702
     DO While ((ITmp.LE.FrozBlock(Iat,0)).AND.(sDihe.LE.0.09d0))
703
        ITmp=ITmp+1
704
        n1=FrozBlock(Iat,Itmp)
705
        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
706
        CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, &
707
             vx4,vy4,vz4,norm4)
708
        val_d=angle_d(vx4,vy4,vz4,norm4,  &
709
             vx5,vy5,vz5,norm5,           &
710
             vx2,vy2,vz2,norm2)
711
        sDihe=abs(sin(val_d))
712
     END DO
713
     IF (debug) WRITE(*,*) 'DBG 4th atom, ITmp, sDihe',ITmp, sDihe
714
     if (sDihe.LE.0.09d0) THEN
715
        !     None of the frozen atoms linked to IAt are good to define the third
716
        !     direction in space...
717
        !     We will look at the other frozen atoms
718
        !     we might improve the search so as to take the atom closest to IAt
719
        ITmp=0
720
        DO I=1,NbAtFrag(IFrag)
721
           JAt=FragAt(IFrag,I)
722
           if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
723
              n1=JAt
724
              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
725
              CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2,   &
726
                   vx4,vy4,vz4,norm4)
727
              val_d=angle_d(vx4,vy4,vz4,norm4,    &
728
                   vx5,vy5,vz5,norm5,              &
729
                   vx2,vy2,vz2,norm2)
730
              if (abs(sin(val_d)).GE.0.09D0) THEN
731
                 ITmp=ITmp+1
732
                 DistFroz(ITmp)=Norm1
733
                 FrozDist(ITmp)=JAt
734
              END IF
735
           END IF
736
        END DO
737
        IF (ITmp.EQ.0) THEN
738
           !     All dihedral angles between frozen atoms are less than 5?
739
           !     (or more than 175?). We have to look at other fragments :-(
740
           DO I=1,NFroz
741
              Jat=Frozen(I)
742
              if (.NOT.DejaFait(JAt)) THEN
743
                 n1=JAt
744
                 CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
745
                 CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2,   &
746
                      vx4,vy4,vz4,norm4)
747
                 val_d=angle_d(vx4,vy4,vz4,norm4,       &
748
                      vx5,vy5,vz5,norm5,                 &
749
                      vx2,vy2,vz2,norm2)
750
                 if (abs(sin(val_d)).GE.0.09D0) THEN
751
                    ITmp=ITmp+1
752
                    DistFroz(ITmp)=Norm1
753
                    FrozDist(ITmp)=JAt
754
                 END IF
755
              END IF
756
           END DO
757
           IF (ITmp.EQ.0) THEN
758
              !     All frozen atoms are in a plane... too bad
759
              WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
760
              WRITE(*,*) 'For now, I do not treat this case'
761
              STOP
762
           END IF
763
        END IF
764
        I1=0
765
        d=1e9
766
        DO I=1,ITmp
767
           IF (DistFroz(I).LE.d) THEN
768
              d=DistFroz(I)
769
              I1=FrozDist(I)
770
           END IF
771
        END DO
772
     ELSE                   !(sDihe.LE.0.09d0)
773
        I1=FrozBlock(IAt,ITmp)
774
     END IF                 !(sDihe.LE.0.09d0)
775
     !     we now have the atom that is closer to the first one and that
776
     !     forms a non 0/Pi dihedral angle
777
     !     ind_zmat(4,1)=I1
778
     !     ind_zmat(4,2)=IAt
779
     !     ind_zmat(4,3)=LiaisonsBis(Iat,1)
780
     !     ind_zmat(4,4)=LiaisonsBis(Iat,2)
781
     n3=LiaisonsBis(Iat,1)
782
     n4=LiaisonsBis(Iat,2)
783
     Call add2indzmat(na,4,I1,Iat,n3,n4,ind_zmat,x,y,z)
784
     LiaisonsBis(Iat,1)=n3
785
     LiaisonsBis(Iat,2)=n4
786
     DejaFait(I1)=.TRUE.
787
     CaFaire(3)=I1
788
     CaFaire(4)=0
789
     IdxCaFaire=4
790
     izm=4
791
     FCaf(I1)=.TRUE.
792

    
793
  CASE(1)
794
     if (debug) WRITE(*,*) 'DBG select case I0 1'
795
     ind_zmat(1,1)=IAt
796
     ind_zmat(2,1)=liaisonsBis(IAt,1)
797
     ind_zmat(2,2)=IAt
798
     DejaFait(IAt)=.TRUE.
799
     DejaFait(LiaisonsBis(Iat,1))=.TRUE.
800
     IdxCaFaire=2
801
     CaFaire(1)=LiaisonsBis(Iat,1)
802
     CaFaire(2)=0
803
     FCaf(LiaisonsBis(Iat,1))=.TRUE.
804

    
805
     !     We search for a third and fourth atoms, first in the FrozBlock atoms
806
     !     It should not be possible to have   (FrozBlock(Iat,0).GT.2) and
807
     !     iat linked to only one atom !       
808

    
809
     IF (FrozBlock(Iat,0).GT.2) THEN
810
        WRITE(*,*) 'I found some inconsistancy: IAt linked to 1'
811
        WRITE(*,*) 'Atom only, but belongs to a frozblock of at '
812
        WRITE(*,*) 'least 3 atoms. IAt, FrozBlock'
813
        WRITE(*,*) Iat,(FrozBlock(IAt,J),J=0,FrozBlock(Iat,0))
814
        STOP
815
     END IF
816

    
817
     !     we calculate the distances between Iat and all other frozen
818
     !     atoms of this fragment, and store only those for which 
819
     !     valence angles are not too close to 0/Pi. (limit:5?)
820

    
821
     ITmp=0
822
     CALL vecteur(LiaisonsBis(Iat,1),IAt,x,y,z,vx2,vy2,vz2,norm2)
823

    
824
     DO I=1,NbAtFrag(IFrag)
825
        JAt=FragAt(IFrag,I)
826
        if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
827
           CALL vecteur(JAt,IAt,x,y,z,vx1,vy1,vz1,norm1)
828
           if (abs(cos(angle(vx1,vy1,vz1,norm1,  &
829
                vx2,vy2,vz2,norm2))).LE.0.996d0) THEN
830
              ITmp=ITmp+1
831
              DistFroz(ITmp)=Norm1
832
              FrozDist(ITmp)=JAt
833
           END IF
834
        END IF
835
     END DO
836

    
837
     IF (ITMP.EQ.0) THEN
838
        !     We have no frozen atoms correct in this fragment, we use
839
        !     atoms from other fragments
840
        DO I=1,NFroz
841
           Jat=Frozen(I)
842
           if (.NOT.DejaFait(JAt)) THEN
843
              CALL vecteur(JAt,IAt,x,y,z,vx1,vy1,vz1,norm1)
844
              if (abs(cos(angle(vx1,vy1,vz1,norm1, &
845
                   vx2,vy2,vz2,norm2))).LE.0.996d0) THEN
846
                 ITmp=ITmp+1
847
                 DistFroz(ITmp)=Norm1
848
                 FrozDist(ITmp)=JAt
849
              END IF
850
           END IF
851
        END DO
852
        IF (ITMP.EQ.0) THEN
853
           WRITE(*,*) 'It seems all frozen atoms are aligned'
854
           WRITE(*,*) 'Case non treated for now :-( '
855
           STOP
856
        END IF
857
     END IF
858

    
859
     I1=0
860
     d=1e9
861
     DO I=1,ITmp
862
        IF (DistFroz(I).LE.d) THEN
863
           I1=FrozDist(I)
864
           d=DistFroz(I)
865
        END IF
866
     END DO
867

    
868
     !     we now have the atom that is closer to the first one and that
869
     !     forms a non 0/Pi valence angle
870
     ind_zmat(3,1)=I1
871
     ind_zmat(3,2)=IAt
872
     ind_zmat(3,3)=LiaisonsBis(Iat,1)
873
     DejaFait(I1)=.TRUE.
874
     CaFaire(2)=I1
875
     FCaf(I1)=.TRUE.
876

    
877

    
878
     !     we search for a fourth atom !
879
     !     We search for a fourth atom, first in the FrozBlock atoms
880
     ITmp=2
881
     sDihe=0.
882
     n2=IAt
883
     n3=LiaisonsBis(Iat,1)
884
     n4=I1
885
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
886
     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
887
     CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2,  &
888
          vx5,vy5,vz5,norm5)
889

    
890
     !     We will look at the other frozen atoms
891
     !     we might improve the search so as to take the atom closest to IAt
892
     ITmp=0
893
     DO I=1,NbAtFrag(IFrag)
894
        JAt=FragAt(IFrag,I)
895
        if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
896
           n1=JAt
897
           CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
898
           CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, &
899
                vx4,vy4,vz4,norm4)
900
           val_d=angle_d(vx4,vy4,vz4,norm4,  &
901
                vx5,vy5,vz5,norm5,   &
902
                vx2,vy2,vz2,norm2)
903
           if (abs(sin(val_d)).GE.0.09D0) THEN
904
              ITmp=ITmp+1
905
              DistFroz(ITmp)=Norm1
906
              FrozDist(ITmp)=JAt
907
           END IF
908
        END IF
909
     END DO
910
     IF (ITmp.EQ.0) THEN
911
        !     All dihedral angles between frozen atoms are less than 5?
912
        !     (or more than 175?). We have to look at other fragments :-(
913
        DO I=1,NFroz
914
           Jat=Frozen(I)
915
           if (.NOT.DejaFait(JAt)) THEN
916
              n1=JAt
917
              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
918
              CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, &
919
                   vx4,vy4,vz4,norm4)
920
              val_d=angle_d(vx4,vy4,vz4,norm4,   &
921
                   vx5,vy5,vz5,norm5,           &
922
                   vx2,vy2,vz2,norm2)
923
              if (abs(sin(val_d)).GE.0.09D0) THEN
924
                 ITmp=ITmp+1
925
                 DistFroz(ITmp)=Norm1
926
                 FrozDist(ITmp)=JAt
927
              END IF
928
           END IF
929
        END DO
930
        IF (ITmp.EQ.0) THEN
931
           !     All frozen atoms are in a plane... too bad
932
           WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
933
           WRITE(*,*) 'For now, I do not treat this case'
934
           STOP
935
        END IF
936
     END IF                 ! ITmp.EQ.0 after scaning fragment
937
     I1=0
938
     d=1e9
939
     DO I=1,ITmp
940
        IF (DistFroz(I).LE.d) THEN
941
           d=DistFroz(I)
942
           I1=FrozDist(I)
943
        END IF
944
     END DO
945

    
946
     !     we now have the atom that is closer to the first one and that
947
     !     forms a non 0/Pi dihedral angle
948
     !     ind_zmat(4,1)=I1
949
     !     ind_zmat(4,2)=IAt
950
     !     ind_zmat(4,3)=ind_zmat(2,1)
951
     !     ind_zmat(4,4)=ind_zmat(3,1)
952
     n3=ind_zmat(2,1)
953
     n4=ind_zmat(3,1)
954
     Call add2indzmat(na,4,I1,IAt,n3,n4,ind_zmat,x,y,z)
955
     ind_zmat(2,1)=n3
956
     ind_zmat(3,1)=n4
957
     DejaFait(I1)=.TRUE.
958
     CaFaire(3)=I1
959
     CaFaire(4)=0
960
     IdxCaFaire=4
961
     izm=4
962
     FCaf(I1)=.TRUE.
963

    
964
  CASE(0)
965
     WRITE(*,*) 'All frozen atoms are separated .. '
966
     WRITE(*,*) 'this case should be treated separately !'
967
     STOP
968
  END SELECT
969

    
970
  if (debug) THEN
971
     WRITE(*,*) 'ind_zmat 1'
972
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
973
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
974
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
975
     DO I=4,izm
976
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2),  &
977
             ind_zmat(I,3), ind_zmat(I,4)
978
     END DO
979
  END IF
980

    
981
  DO I=1,izm
982
     Idx_zmat(ind_zmat(I,1))=i
983
  END Do
984

    
985
  !     at least first four (frozen) atoms of this fragment done...
986
  !     we empty the 'cafaire' array before pursuing
987
  IAFaire=1
988
  DO WHILE (CaFaire(IaFaire).NE.0)
989
     n1=CaFaire(IaFaire)
990
     n2=ind_zmat(idx_zmat(N1),2)
991
     if (idx_zmat(N1).EQ.2) THEN
992
        !     We have a (small) problem: we have to add atoms linked to the 2nd
993
        !     atom of the zmat. This is a pb because we do not know
994
        !     which atom to use to define the dihedral angle
995
        !     we take the third atom of the zmat
996
        n3=ind_zmat(3,1)
997
     ELSE
998
        n3=ind_zmat(idx_zmat(n1),3)
999
     END IF
1000
     IF (LiaisonsBis(n1,0).GE.1) THEN
1001
        IAt=LiaisonsBis(n1,1)
1002
        if (.NOT.DejaFait(Iat)) THEN
1003
           izm=izm+1
1004
           if (debug) WRITE(*,*) ">0< Adding atom ",Iat,"position izm=",izm
1005
           !           ind_zmat(izm,1)=iat
1006
           !           ind_zmat(izm,2)=n1
1007
           !           ind_zmat(izm,3)=n2
1008
           !           ind_zmat(izm,4)=n3
1009
           Call add2indzmat(na,izm,iat,n1,n2,n3,ind_zmat,x,y,z)
1010
           Idx_zmat(iat)=izm
1011
           n3=iat
1012
           IF (.NOT.FCaf(Iat)) THEN
1013
              CaFaire(IdxCaFaire)=iat
1014
              IdxCaFaire=IdxCaFaire+1
1015
              CaFaire(IdxCaFaire)=0
1016
              FCaf(Iat)=.TRUE.
1017
           END IF
1018
           DejaFait(iat)=.TRUE.
1019
        END IF
1020
     END IF
1021
     DO I=2,LiaisonsBis(n1,0)
1022
        IAt=LiaisonsBis(n1,I)
1023
! PFL 29.Aug.2008 ->
1024
! We dissociate the test on 'DejaFait' that indicates that this atom
1025
! has already been put in the Zmat
1026
! from the test on FCaf that indicates that this atom has been put in the
1027
! 'CAFaire' list that deals with identifying its connectivity.
1028
! Those two test might differ for a frozen atom linked to non frozen atoms.
1029
        IF (.NOT.DejaFait(Iat)) THEN
1030
           izm=izm+1
1031
           if (debug) WRITE(*,*) ">1< Adding atom ",Iat,"position izm=",izm
1032
           !           ind_zmat(izm,1)=iat
1033
           !           ind_zmat(izm,2)=n1
1034
           !           ind_zmat(izm,3)=n2
1035
           !           ind_zmat(izm,4)=n3
1036
           Call add2indzmat(na,izm,iat,n1,n2,n3,ind_zmat,x,y,z)
1037
           idx_zmat(iat)=izm
1038
           DejaFait(iat)=.TRUE.
1039
        END IF
1040
        IF (.NOT.FCaf(Iat)) THEN
1041
           CaFaire(IdxCaFaire)=iat
1042
           IdxCaFaire=IdxCaFaire+1
1043
           CaFaire(IdxCaFaire)=0
1044
           FCaf(Iat)=.TRUE.
1045
        END IF
1046
! <- PFL 29.Aug.2008
1047
     END DO
1048
     IaFaire=IaFaire+1
1049
  END Do                    ! DO WHILE CaFaire
1050

    
1051
  if (debug) THEN
1052
     WRITE(*,*) 'ind_zmat 2, izm=',izm
1053
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
1054
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
1055
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
1056
     DO I=4,izm
1057
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2),  &
1058
             ind_zmat(I,3), ind_zmat(I,4)
1059
     END DO
1060
  END IF
1061

    
1062
  !     We have finished putting atoms linked to the first one
1063
  !     we will add other frozen atoms of this fragment
1064
  DO I=1,NbAtFrag(IFrag)
1065
     Iat=FragAt(IFrag,I)
1066
     if (debug) WRITE(*,*) "DBG: I,Iat,Frozat,dejafait,frozbloc",I,Iat,FrozAt(Iat), DejaFait(Iat),FrozBlock(Iat,0)
1067
     IF (FrozAt(Iat).AND.(.NOT.DejaFait(Iat))) THEN
1068
        !     in order to have the zmat as connected as possible,
1069
        !     we look if this atom belongs to a frozblock
1070
        if (debug) WRITE(*,*) 'Treating atom I,Iat, FrozBlock ',I,Iat, FrozBlock(Iat,0)
1071
        IF (FrozBlock(Iat,0).EQ.1) THEN
1072
           !     it is a lonely atom :-)
1073
           d=1e9
1074
           DO J=1,izm
1075
              Call vecteur(iat,ind_zmat(j,1),x,y,z,vx1,vy1,vz1,norm1)
1076
              if (norm1.le.d) THEN
1077
                 Jat=j
1078
                 d=norm1
1079
              END IF
1080
           END DO
1081
           n2=ind_zmat(jat,1)
1082
           SELECT CASE(jat)
1083
           CASE (1)
1084
              n3=ind_zmat(2,1)
1085
              n4=ind_zmat(3,1)
1086
           CASE (2)
1087
              n3=ind_zmat(1,1)
1088
              n4=ind_zmat(3,1)
1089
           CASE DEFAULT
1090
              n3=ind_zmat(jAt,2)
1091
              n4=ind_zmat(jat,3)
1092
           END SELECT
1093
           izm=izm+1
1094
           !           ind_zmat(izm,1)=iat
1095
           !           ind_zmat(izm,2)=n2
1096
           !           ind_zmat(izm,3)=n3
1097
           !           ind_zmat(izm,4)=n4
1098
           Call add2indzmat(na,izm,iat,n2,n3,n4,ind_zmat,x,y,z)
1099
           DejaFait(iat)=.TRUE.
1100
           idx_zmat(iat)=iat
1101
        ELSE        !(FrozBlock(Iat,0).EQ.1)
1102
           !     we have a group of atoms
1103
           !     We search for the atom of the group with the most links
1104
           ITmp=-1
1105
           DO J=1,FrozBlock(Iat,0)
1106
              Jat=FrozBlock(Iat,J)
1107
              IF ((.NOT.DejaFait(Jat)).AND.  &
1108
                   (LiaisonsBis(Jat,0).GT.ITMP)) THEN
1109
                 JL=Jat
1110
                 ITmp=LiaisonsBis(Jat,0)
1111
              END IF
1112
           END DO
1113
           if (debug) WRITE(*,*) 'JL,ITmp:',JL,ITmp
1114
           Iat=JL
1115
           d=1e9
1116
           DO J=1,izm
1117
              Call vecteur(iat,ind_zmat(j,1),x,y,z, vx1,vy1,vz1,norm1)
1118
              if (norm1.le.d) THEN
1119
                 Jat=j
1120
                 d=norm1
1121
              END IF
1122
           END DO
1123
           if (debug) WRITE(*,*) 'Jat,d:',Jat,d
1124
           n2=ind_zmat(jat,1)
1125
           SELECT CASE(jat)
1126
           CASE (1)
1127
              n3=ind_zmat(2,1)
1128
              n4=ind_zmat(3,1)
1129
           CASE (2)
1130
              n3=ind_zmat(1,1)
1131
              n4=ind_zmat(3,1)
1132
           CASE DEFAULT
1133
              n3=ind_zmat(jAt,2)
1134
              n4=ind_zmat(jat,3)
1135
           END SELECT
1136
           izm=izm+1
1137
           !           ind_zmat(izm,1)=iat
1138
           !           ind_zmat(izm,2)=n2
1139
           !           ind_zmat(izm,3)=n3
1140
           !           ind_zmat(izm,4)=n4
1141
           Call add2indzmat(na,izm,iat,n2,n3,n4,ind_zmat,x,y,z)
1142
           idx_zmat(iat)=izm
1143
           DejaFait(iat)=.TRUE.
1144
           IdxCaFaire=2
1145
           CaFaire(1)=iat
1146
           CaFaire(2)=0
1147
           FCaf(Iat)=.TRUE.
1148

    
1149
           if (debug) WRITE(*,*) izm,iat,n2,n3,n4
1150

    
1151
           IaFaire=1
1152
           DO WHILE (CaFaire(IaFaire).NE.0)
1153
              n1=CaFaire(IaFaire)
1154
              n2=ind_zmat(idx_zmat(N1),2)
1155
              if (debug) WRITE(*,*) 'DBG Cafaire, Iafaire,n1,n2',Iafaire,n1,n2
1156
              if (idx_zmat(N1).EQ.2) THEN
1157
                 !     We have a (small) problem: we have to add atoms linked to the 2nd
1158
                 !     atom of the zmat. This is a pb because we do not know
1159
                 !     which atom to use to define the dihedral angle
1160
                 !     we take the third atom of the zmat
1161
                 n3=ind_zmat(3,1)
1162
              ELSE
1163
                 n3=ind_zmat(idx_zmat(n1),3)
1164
              END IF
1165
              if (debug) WRITe(*,*) 'DBG :n3,liaisonBis',n3,LiaisonsBis(n1,0)
1166
              DO I3=1,LiaisonsBis(n1,0)
1167
! PFL 29.Aug.2008 ->
1168
! We dissociate the test on 'DejaFait' that indicates that this atom
1169
! has already been put in the Zmat
1170
! from the test on FCaf that indicates that this atom has been put in the
1171
! 'CAFaire' list that deals with identifying its connectivity.
1172
! Those two test might differ for a frozen atom linked to non frozen atoms.
1173
                 IAt=LiaisonsBis(n1,I3)
1174
                 if (.NOT.DejaFait(Iat)) THEN
1175
                    izm=izm+1
1176
                    !                 ind_zmat(izm,1)=iat
1177
                    !                 ind_zmat(izm,2)=n1
1178
                    !                 ind_zmat(izm,3)=n2
1179
                    !                 ind_zmat(izm,4)=n3
1180
                    Call add2indzmat(na,izm,iat,n1,n2,n3,ind_zmat,x,y,z)
1181
                    idx_zmat(iat)=izm
1182
                    if (I3.EQ.1) n3=ind_zmat(izm,1)
1183
                    DejaFait(Iat)=.TRUE.
1184
                 END IF
1185
                 If (.NOT.FCaf(Iat)) THEN
1186
                    CaFaire(IdxCaFaire)=iat
1187
                    IdxCaFaire=IdxCaFaire+1
1188
                    CaFaire(IdxCaFaire)=0
1189
                    FCaf(Iat)=.TRUE.
1190
                 END IF
1191
! <- PFL 29.Aug.2008
1192
              END DO
1193
              IaFaire=IaFaire+1
1194
           END Do           ! DO WHILE CaFaire
1195
        END IF
1196
     END IF
1197
     if (debug) THEN
1198
        WRITE(*,*) 'ind_zmat 3'
1199
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
1200
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
1201
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
1202
        DO Ip=4,izm
1203
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1), ind_zmat(Ip,2),  &
1204
                ind_zmat(Ip,3), ind_zmat(Ip,4)
1205
        END DO
1206
     END IF
1207

    
1208
  END DO
1209

    
1210
  FrozFrag(IFrag,1)=-1
1211
  FrozFrag(IFrag,2)=-1
1212
  FrozFrag(IFrag,3)=-1
1213

    
1214
  !     we start again with the rest of the molecule...
1215
  ! v 1.01 We add the fragment in the order corresponding to FrozFrag
1216
  KMax=0
1217
  DO I=1,NbFrag
1218
     IF (FrozFrag(I,1).NE.0) KMax=KMax+1
1219
  END DO
1220

    
1221
  IF (DEBUG) WRITE(*,*) "Adding the",Kmax,"fragments with frozen atoms"
1222
  DO K=1, KMax-1
1223
     IFrag=0
1224
     I0=0
1225
     I1=0
1226
     DO I=1,NbFrag
1227
        SELECT CASE(FrozFrag(I,3)-I0)
1228
        CASE (1:)
1229
           IFrag=I
1230
           I0=FrozFrag(I,3)
1231
           I1=FrozFrag(I,2)
1232
        CASE (0)
1233
           if (FrozFrag(I,2).GT.I1) THEN
1234
              IFrag=I
1235
              I0=FrozFrag(I,3)
1236
              I1=FrozFrag(I,2)
1237
           END IF
1238
        END SELECT
1239
     END DO
1240

    
1241
     FrozFrag(IFrag,1)=-1
1242
     FrozFrag(IFrag,2)=-1
1243
     FrozFrag(IFrag,3)=-1
1244

    
1245
     if (debug) WRITE(*,*) "Adding Fragment ",IFrag,K
1246

    
1247
     DO I=1,NbAtFrag(IFrag)
1248
        Iat=FragAt(IFrag,I)
1249
        IF (FrozAt(Iat).AND.(.NOT.DejaFait(Iat))) THEN
1250
           !     in order to have the zmat as connected as possible,
1251
           !     we look if this atom belongs to a frozblock
1252
           IF (FrozBlock(Iat,0).EQ.1) THEN
1253
              !     it is a lonely atom :-)
1254
              if (debug) write(*,*) "DBG 4: Lonely atom",Iat
1255
              d=1e9
1256
              DO J=1,izm
1257
                 Call vecteur(iat,ind_zmat(j,1),x,y,z,vx1,vy1,vz1,norm1)
1258
                 if (norm1.le.d) THEN
1259
                    Jat=j
1260
                    d=norm1
1261
                 END IF
1262
              END DO
1263
              n2=ind_zmat(jat,1)
1264
              SELECT CASE(jat)
1265
              CASE (1)
1266
                 n3=ind_zmat(2,1)
1267
                 n4=ind_zmat(3,1)
1268
              CASE (2)
1269
                 n3=ind_zmat(1,1)
1270
                 n4=ind_zmat(3,1)
1271
              CASE DEFAULT
1272
                 n3=ind_zmat(jAt,2)
1273
                 n4=ind_zmat(jat,3)
1274
              END SELECT
1275
              izm=izm+1
1276
              !              ind_zmat(izm,1)=iat
1277
              !              ind_zmat(izm,2)=n2
1278
              !              ind_zmat(izm,3)=n3
1279
              !              ind_zmat(izm,4)=n4
1280
              Call add2indzmat(na,izm,iat,n2,n3,n4,ind_zmat,x,y,z)
1281
              idx_zmat(iat)=izm
1282
              DejaFait(iat)=.TRUE.
1283
           ELSE
1284
              !     we have a group of atoms
1285
              !     We search for the atom of the group with the most links
1286
              if (debug) write(*,*) "DBG 4b: ",Iat," belong to frozblock",FrozBlock(Iat,0)
1287
              ITmp=-1
1288
              DO J=1,FrozBlock(Iat,0)
1289
                 Jat=FrozBlock(Iat,J)
1290
                 IF ((.NOT.DejaFait(Jat)).AND.         &
1291
                      (LiaisonsBis(Jat,0).GT.ITMP)) THEN
1292
                    JL=Jat
1293
                    ITmp=LiaisonsBis(Jat,0)
1294
                 END IF
1295
              END DO
1296
              Iat=JL
1297
              d=1e9
1298
              DO J=1,izm
1299
                 Call vecteur(iat,ind_zmat(j,1),x,y,z,vx1,vy1,vz1,norm1)
1300
                 if (norm1.le.d) THEN
1301
                    Jat=j
1302
                    d=norm1
1303
                 END IF
1304
              END DO
1305
              n2=ind_zmat(jat,1)
1306
              SELECT CASE(jat)
1307
              CASE (1)
1308
                 n3=ind_zmat(2,1)
1309
                 n4=ind_zmat(3,1)
1310
              CASE (2)
1311
                 n3=ind_zmat(1,1)
1312
                 n4=ind_zmat(3,1)
1313
              CASE DEFAULT
1314
                 n3=ind_zmat(jAt,2)
1315
                 n4=ind_zmat(jat,3)
1316
              END SELECT
1317
              izm=izm+1
1318
              !              ind_zmat(izm,1)=iat
1319
              !              ind_zmat(izm,2)=n2
1320
              !              ind_zmat(izm,3)=n3
1321
              !              ind_zmat(izm,4)=n4
1322
              Call add2indzmat(na,izm,Iat,n2,n3,n4,ind_zmat,x,y,z)
1323
              idx_zmat(iat)=izm
1324
              DejaFait(iat)=.TRUE.
1325
              IdxCaFaire=2
1326
              CaFaire(1)=iat
1327
              CaFaire(2)=0
1328
              FCaf(Iat)=.TRUE.
1329
              IaFaire=1
1330
              DO WHILE (CaFaire(IaFaire).NE.0)
1331
                 n1=CaFaire(IaFaire)
1332
                 n2=ind_zmat(idx_zmat(N1),2)
1333
                 if (idx_zmat(N1).EQ.2) THEN
1334
                    !     We have a (small) problem: we have to add atoms linked to the 2nd
1335
                    !     atom of the zmat. This is a pb because we do not know
1336
                    !     which atom to use to define the dihedral angle
1337
                    !     we take the third atom of the zmat
1338
                    n3=ind_zmat(3,1)
1339
                 ELSE
1340
                    n3=ind_zmat(idx_zmat(n1),3)
1341
                 END IF
1342
                 DO I3=1,LiaisonsBis(n1,0)
1343
                    IAt=LiaisonsBis(n1,I3)
1344
! PFL 29.Aug.2008 ->
1345
! We dissociate the test on 'DejaFait' that indicates that this atom
1346
! has already been put in the Zmat
1347
! from the test on FCaf that indicates that this atom has been put in the
1348
! 'CAFaire' list that deals with identifying its connectivity.
1349
! Those two test might differ for a frozen atom linked to non frozen atoms.
1350
                    IF (.NOT.DejaFait(Iat)) THEN
1351
                       izm=izm+1
1352
                       Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
1353
                       idx_zmat(iat)=izm
1354
                       n3=iat
1355
                       DejaFait(Iat)=.TRUE.
1356
                    END IF
1357
                    IF (.NOT.FCaf(Iat)) THEN
1358
                       CaFaire(IdxCaFaire)=iat
1359
                       IdxCaFaire=IdxCaFaire+1
1360
                       CaFaire(IdxCaFaire)=0
1361
                       FCaf(Iat)=.TRUE.
1362
                    END IF
1363
! <- PFL 29.Aug.2008
1364
                 END DO
1365
                 IaFaire=IaFaire+1
1366
              END Do              ! DO WHILE CaFaire
1367
           END IF
1368
        END IF
1369

    
1370
        if (debug) THEN
1371
           WRITE(*,*) 'ind_zmat 4'
1372
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
1373
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
1374
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
1375
           DO Ip=4,izm
1376
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1), ind_zmat(Ip,2), &
1377
                   ind_zmat(Ip,3), ind_zmat(Ip,4)
1378
           END DO
1379
        END IF
1380

    
1381
     END DO ! loop on atoms of fragment IFrag
1382
  END DO                    ! Loop on all fragments
1383

    
1384
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1385
!
1386
!        General case
1387
!
1388
! 2 - we have all frozen atoms done... now comes the non frozen atoms
1389
! and we should fragment the fragments !
1390
!
1391
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1392

    
1393
! First, we get rid of bonds between frozen atoms
1394
! the trick is to do this on liaisons while keeping LiaisonsBis ok.
1395

    
1396
  if (debug) THEN
1397
     WRITE(*,*) 'Frozen',NFroz
1398
     WRITE(*,'(20(1X,I3))') (Frozen(I),I=1,NFroz)
1399
  END IF
1400

    
1401
  DO I=1,na
1402
     DO J=0,NMAxL
1403
        LiaisonsIni(I,J)=LiaisonsBis(I,J)
1404
     END DO
1405
! PFL 29.Aug.2008 -> 
1406
! We re-initialize FCaf in order to add frozen atoms that are connected
1407
! to non frozen atoms
1408
     FCaf(I)=.FALSE.
1409
! <- 29.Aug.2008
1410
  END DO
1411

    
1412
  DO I=1,Na
1413
     IF (FrozAt(I)) THEN
1414
        Iat=I
1415
        if (debug) WRITE(*,'(20(1X,I3))') I,Iat,  &
1416
             (LiaisonsIni(Iat,Itmp),ITmp=0,LiaisonsIni(Iat,0))
1417
        DO J=1,LiaisonsIni(Iat,0)
1418
           Jat=LiaisonsIni(Iat,J)
1419
           Call Annul(Liaisons,Iat,Jat)
1420
           Call Annul(Liaisons,Jat,Iat)
1421
           Call Annul(LiaisonsIni,Jat,Iat)
1422
           Liaisons(Iat,0)=Liaisons(Iat,0)-1
1423
           Liaisons(Jat,0)=Liaisons(Jat,0)-1
1424
           LiaisonsIni(Jat,0)=LiaisonsIni(Jat,0)-1
1425
        END DO
1426
     END IF
1427
  END DO
1428

    
1429
  if (debug) THEN
1430
     WRITE(*,*) "Connections calculees"
1431
     DO IL=1,na
1432
        WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
1433
     END DO
1434
  END IF
1435

    
1436

    
1437

    
1438
  ! We re-decompose the system into fragments, but we omit on purpuse
1439
  ! fragments consisting only of frozen atoms
1440
  ! Now, frozblock will contain the frozen atom indices in a given fragment
1441

    
1442
  DO I=1,na
1443
     Fragment(I)=0 
1444
     NbAtFrag(I)=0
1445
     FrozBlock(I,0)=0
1446
  END DO
1447
  IdxFrag=0
1448
  NbFrag=0
1449

    
1450
  DO I=1,na
1451
     IdxCAfaire=1
1452
     CaFaire(IdxCaFaire)=0
1453

    
1454
     if (debug) WRITE(*,*) 'Treating atom I, fragment(I)',I,Fragment(I)
1455
     IF (.NOT.FrozAt(I).OR.(Liaisons(I,0).NE.0)) THEN 
1456
        IF (Fragment(I).EQ.0) THEN
1457
           IdxFrag=IdxFrag+1
1458
           NbFrag=NbFrag+1
1459
           IFrag=IdxFrag
1460
           Fragment(I)=IFrag
1461
           NbAtFrag(IFrag)=NbAtFrag(IFrag)+1
1462
           FragAt(IFrag,NbAtFrag(IFrag))=I
1463
        ELSE  ! fragment(I).EQ.0
1464
           IFrag=Fragment(I)
1465
        END IF     ! fragment(I).EQ.0
1466
        DO J=1,Liaisons(I,0)
1467
           Iat=Liaisons(I,J)
1468
           IF ((Fragment(Iat).NE.0).AND.(Fragment(Iat).NE.IFrag)) THEN
1469
              WRITE(*,*) 'Error : Atoms ',I,' and ',Iat
1470
              WRITE(*,*) 'are linked, but belongs to fragments ',IFrag,' and ',Fragment(Iat)
1471
              STOP
1472
           END IF
1473
           IF (Fragment(Iat).EQ.0) THEN
1474
              IF (.NOT.FCaf(IAt)) THEN
1475
                 CaFaire(IdxCaFaire)=Iat
1476
                 IdxCAfaire=IdxCAFaire+1
1477
                 FCaf(IAt)=.TRUE.
1478
              END IF
1479
              Fragment(Iat)=IFrag
1480
              NbAtFrag(IFrag)=NbAtFrag(IFrag)+1
1481
              FragAt(IFrag,NbAtFrag(IFrag))=Iat
1482
           END IF
1483
        END DO
1484
        CaFaire(IdxCaFaire)=0
1485

    
1486
        If (debug) WRITE(*,'(1X,A,1000I4)') 'Cafaire:',CaFaire(1:IdxCaFaire)
1487
        If (debug) WRITE(*,*) 'IFrag=',IFrag
1488

    
1489
        IAfaire=1
1490
        DO WHILE (CaFaire(IAfaire).NE.0)
1491
           Iat=CaFaire(IaFaire)
1492
           IF (.NOT.FrozAt(Iat).OR.(Liaisons(Iat,0).NE.0)) THEN 
1493
              if (debug) WRITE(*,*) 'Cafaire treating ',Iat
1494
              IF (Fragment(Iat).EQ.0) THEN
1495
                 WRITE(*,*) 'Error: Atom ',Iat,' does not belong to any fragment !'
1496
                 STOP
1497
              END IF
1498

    
1499
              DO J=1,Liaisons(Iat,0)
1500
                 Jat=Liaisons(Iat,J)
1501
                 IF ((Fragment(Jat).NE.0).AND.(Fragment(Jat).NE.IFrag)) THEN
1502
                    WRITE(*,*) 'Error: Atoms ',Iat,' and ',Liaisons(Iat,J)
1503
                    WRITE(*,*) 'are linked, but belongs to fragments ',IFrag,' and ',Fragment(Liaisons(Iat,J))
1504
                    STOP
1505
                 END IF
1506
                 IF (Fragment(Jat).EQ.0) THEN
1507
                    IF (.NOT.FCaf(Liaisons(Iat,J))) THEN
1508
                       CaFaire(IdxCaFaire)=Liaisons(Iat,J)
1509
                       IdxCAfaire=IdxCAFaire+1
1510
                       FCaf(Liaisons(Iat,J))=.TRUE.
1511
                    END IF
1512
                    Fragment(Jat)=IFrag
1513
                    NbAtFrag(IFrag)=NbAtFrag(IFrag)+1
1514
                    FragAt(IFrag,NbAtFrag(IFrag))=Jat
1515
                 END IF
1516
              END DO
1517
              CaFaire(IdxCaFaire)=0
1518
              If (debug) WRITE(*,*) 'IAfaire, IdxCaFaire,Cafaire :',IAfaire,IdxCafaire,' == ',CaFaire(IaFaire+1:IdxCaFaire)
1519
              IAFaire=IAFaire+1
1520
           END IF
1521
        END DO
1522
     END IF
1523
  END DO
1524

    
1525
  ! We compute FrozBlock now
1526
  DO IFrag=1,NbFrag
1527
     DO I=1,NbAtFrag(IFrag)
1528
        Iat=FragAt(IFrag,I)
1529
        IF (FrozAt(Iat)) THEN
1530
           FrozBlock(IFrag,0)=FrozBlock(IFrag,0)+1
1531
           FrozBlock(IFrag,FrozBlock(IFrag,0))=IAt
1532
        END IF
1533
     END DO
1534
  END DO
1535

    
1536

    
1537
  if (debug) THEN
1538
     WRITE(*,*) 'I found ',NbFrag, 'fragments with respectively '
1539
     WRITE(*,*) (NbAtFrag(I),I=1,NbFrag), 'atoms'
1540
     WRITE(*,*) "Fragments atoms are now listed as F in the following"
1541
     DO I=1,NbFrag
1542
        WRITE(*,*) Na
1543
        WRITE(*,*) 'Fragment ', i
1544
        DO J=1,Na
1545
           AtName=Nom(Atome(J))
1546
           IF (Fragment(J).EQ.I) AtName='F'
1547
           WRITE(*,'(1X,A5,3(1X,F10.6))') AtName,X(J),Y(J),Z(J)
1548
        END DO
1549
        !     WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
1550
     END DO
1551

    
1552
     DO I=1,NbFrag
1553
        WRITE(*,*) 'Fragment ', i
1554
        WRITE(*,'(A,20(1X,I3))') "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
1555
        WRITE(*,'(A,20(1X,I3))') "FrozBlock:",(FrozBlock(I,j),j=0,FrozBlock(I,0))
1556
     END DO
1557
  END IF
1558

    
1559

    
1560
  NFroz=0
1561
  DO IFrag=1,NbFrag
1562
     If (FrozBlock(IFrag,0).NE.0) NFroz=NFroz+1
1563
  END DO
1564

    
1565
  IF (DEBUG) WRITE(*,'(1X,A,I4,A,I4,A)') "Found ",NFroz," fragmenst with frozen atoms, out of",NbFrag," fragments"
1566

    
1567
  ! We will now add the fragments containing non frozen atoms.
1568
  ! I am not sure that there is a best strategy to add them...
1569
  ! so we add them in the order they were created :-(
1570
  ! First only block with frozen atoms are added
1571
  izm0=-1
1572
  IFrag=1
1573
  FCaf=.FALSE.
1574

    
1575
  DO IFragFroz=1,NFroz
1576
     DO WHILE (FrozBlock(IFrag,0).EQ.0)
1577
        IFrag=IFrag+1
1578
     END DO
1579
     ! each atom has at least one frozen atom that will serve as the seed
1580
     ! of this fragment.
1581
     if (debug) WRITE(*,*) 'Adding fragment :', ifrag,'with',FrozBlock(IFrag,0),' frozen atoms'
1582
     IF (FrozBlock(IFrag,0).EQ.1) THEN
1583
        ! There is only one frozen atom, easy case ...
1584
        Iat=FrozBlock(IFrag,1)
1585
        if (debug) WRITE(*,*) 'Only frozen atom of frag',iat
1586
        DejaFait(iat)=.TRUE.
1587
        IdxCaFaire=2
1588
        CaFaire(1)=iat
1589
        CaFaire(2)=0
1590
        FCaf(Iat)=.TRUE.
1591
        IaFaire=1
1592
        DO WHILE (CaFaire(IaFaire).NE.0)
1593
           n1=CaFaire(IaFaire)
1594
           SELECT CASE(idx_zmat(n1))
1595
           CASE (1)
1596
              n2=ind_zmat(2,1)
1597
              n3=ind_zmat(3,1)
1598
           CASE (2)
1599
              n2=ind_zmat(1,1)
1600
              n3=ind_zmat(3,1)
1601
           CASE (3:)
1602
              n2=ind_zmat(idx_zmat(n1),2)
1603
              n3=ind_zmat(idx_zmat(n1),3)
1604
           END SELECT
1605

    
1606
           DO I3=1,Liaisons(n1,0)
1607
              IAt=Liaisons(n1,I3)
1608
              ! PFL 2007.Apr.16 ->
1609
              ! We add ALL atoms linked to n1 to CaFaire
1610
              ! But we delete their link to n1
1611
              IF (.NOT.FCaf(Iat)) THEN
1612
                 CaFaire(IdxCaFaire)=Iat
1613
                 IdxCaFaire=IdxCaFaire+1
1614
                 CaFaire(IdxCaFaire)=0
1615
              END IF
1616
              Call Annul(Liaisons,Iat,n1)
1617
              Liaisons(iat,0)=Liaisons(Iat,0)-1
1618
              ! <- PFL 2007.Apr.16
1619
              IF (.NOT.DejaFait(Iat)) THEN
1620
                 ! we add it to the Zmat
1621
                 izm=izm+1
1622
                 !              ind_zmat(izm,1)=iat
1623
                 !              ind_zmat(izm,2)=n1
1624
                 !              ind_zmat(izm,3)=n2
1625
                 !              ind_zmat(izm,4)=n3
1626
                 Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
1627
                 idx_zmat(iat)=izm
1628
                 !                  Call Annul(Liaisons,n1,iat)
1629
                 n3=iat
1630
                 DejaFait(Iat)=.TRUE.
1631
              END IF
1632
           END DO
1633
              IaFaire=IaFaire+1
1634

    
1635
              if (debug) THEN
1636
                 WRITE(*,*) 'ind_zmat 5'
1637
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
1638
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
1639
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
1640
                 DO Ip=4,izm
1641
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
1642
                         ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
1643
                 END DO
1644
              END IF
1645
              
1646
           END Do              ! DO WHILE CaFaire
1647

    
1648

    
1649
        ELSE     ! FrozBlock(I,0)==1
1650
           if (debug) WRITE(*,*) 'IFrag=',IFrag,'Frozblock(I,0)/=1',Frozblock(IFrag,0)
1651
           ! We have many frozen atoms for one fragment. We will have to choose
1652
           ! the first one, and then to decide when to swich...
1653
           Iat=0
1654
           I0=-1
1655
           DO I=1,FrozBlock(IFrag,0)
1656
              JAt=FrozBlock(IFrag,I)
1657
              if (debug) WRITE(*,*) Jat, &
1658
                   (LiaisonsBis(Jat,Itmp),Itmp=0,LiaisonsBis(Jat,0))
1659
              IF (LiaisonsBis(Jat,0).GT.I0) THEN
1660
                 I0=LiaisonsBis(Jat,0)
1661
                 Iat=Jat
1662
              END IF
1663
           END DO
1664
           ! Iat is the starting frozen atom
1665
           IF (debug) WRITE(*,*) 'Choosing atom ',iat,'as a starting atom for frag',ifrag
1666
           DejaFait(iat)=.TRUE.
1667
           IdxCaFaire=2
1668
           CaFaire(1)=iat
1669
           CaFaire(2)=0
1670
           IaFaire=1
1671
           FCaf(Iat)=.TRUE.
1672
           DO WHILE (CaFaire(IaFaire).NE.0)
1673
              n1=CaFaire(IaFaire)
1674
              if (debug) WRITE(*,*) 'Iafaire,n1,dejafait,idx_zmat', &
1675
                   IaFaire,n1,    DejaFait(n1),idx_zmat(n1)
1676
              if (debug) WRITE(*,*) 'Cafaire', &
1677
                   (CaFaire(J),J=Iafaire,IdxCAfaire)
1678
              SELECT CASE(idx_zmat(n1))
1679
              CASE (1)
1680
                 n2=ind_zmat(2,1)
1681
                 n3=ind_zmat(3,1)
1682
              CASE (2)
1683
                 n2=ind_zmat(1,1)
1684
                 n3=ind_zmat(3,1)
1685
              CASE (3:)
1686
                 n2=ind_zmat(idx_zmat(n1),2)
1687
                 n3=ind_zmat(idx_zmat(n1),3)
1688
              END SELECT
1689

    
1690
              if (debug) WRITE(*,*) "DBG n1,Liaisons(n1,0)",n1,Liaisons(n1,0)
1691
              DO I3=1,Liaisons(n1,0)
1692
                 IAt=Liaisons(n1,I3)                 
1693
                 if (debug) WRITE(*,*) "DBG I3,Iat",I3,Iat
1694
                 ! PFL 2007.Apr.16 ->
1695
                 ! We add ALL atoms linked to n1 to CaFaire
1696
                 ! But we delete their link to n1
1697
!! PFL 2007.Aug.28 ->
1698
! re-add the test on FCaf in order not to put the same atom more than once in 
1699
! CaFaire array.
1700
                 IF (.NOT.FCaf(Iat)) THEN
1701
                    CaFaire(IdxCaFaire)=Iat
1702
                    IdxCaFaire=IdxCaFaire+1
1703
                    CaFaire(IdxCaFaire)=0
1704
                    FCaf(Iat)=.TRUE.
1705
                    if (debug) WRITE(*,*) 'Adding ',iat,' to CaFaire. I3=',I3,'IdxCafaire=',IdxCafaire
1706
                    if (debug) WRITE(*,*) 'Cafaire',(CaFaire(J),J=Iafaire,IdxCAfaire)
1707

    
1708
                 END IF
1709
!! <-- PFL 2007.Aug.28
1710

    
1711
                 Call Annul(Liaisons,Iat,n1)
1712
                 Liaisons(iat,0)=Liaisons(Iat,0)-1
1713
                 ! <- PFL 2007.Apr.16
1714
                 IF (.NOT.DejaFait(Iat)) THEN
1715
                    izm=izm+1
1716
                    !                 ind_zmat(izm,1)=iat
1717
                    !                 ind_zmat(izm,2)=n1
1718
                    !                 ind_zmat(izm,3)=n2
1719
                    !                 ind_zmat(izm,4)=n3
1720
                    Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
1721
                    idx_zmat(iat)=izm
1722
                    !                  Call Annul(Liaisons,n1,iat)
1723

    
1724
                    n3=iat
1725
                    DejaFait(Iat)=.TRUE.
1726
                 END IF
1727
              END DO
1728

    
1729
              if (debug) WRITE(*,*) 'I3,n1,Liaisons(n1,:)', &
1730
                   I3,n1,(Liaisons(n1,J),J=0,Liaisons(n1,0))
1731

    
1732

    
1733
              if (debug.AND.(izm.GT.izm0)) THEN
1734
                 WRITE(*,*) 'ind_zmat 6, izm=',izm
1735
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
1736
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
1737
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2),  &
1738
                      ind_zmat(3,3)
1739
                 DO Ip=4,izm
1740
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1), &
1741
                         ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
1742
                 END DO
1743
                 izm0=izm
1744
              END IF
1745

    
1746
              IaFaire=IaFaire+1
1747

    
1748

    
1749
           END Do              ! DO WHILE CaFaire
1750

    
1751
        END IF  ! FrozBlock(I,0)==1
1752

    
1753
        IFrag=IFrag+1
1754
     END DO                    ! Loop on all fragments
1755

    
1756

    
1757
     DO IFrag=1,NbFrag
1758
        IF (FrozBlock(IFrag,0).EQ.0) THEN
1759
           if (debug) WRITE(*,*) 'Adding fragment :', ifrag,'with no frozen atoms'
1760
           ! We have no frozen atoms for this fragment. We will have to choose
1761
           ! the first atom to put.
1762
           ! For now, we do not care of the 'central' atom (ie the one with
1763
           ! no CP atoms...)
1764
           ! We just take the atom that is the closest to the already placed
1765
           ! atoms !
1766

    
1767

    
1768
           I0=0
1769
           I1=0
1770
           d=1e9
1771
           DO J=1,izm
1772
              Jat=ind_zmat(J,1)
1773
              DO I=1,NbAtfrag(IFrag)
1774
                 IAt=FragAt(IFrag,I)
1775
                 IF (.NOT.DejaFait(Iat)) THEN
1776
                    Call vecteur(Iat,Jat,x,y,z,vx1,vy1,vz1,norm1)
1777
                    IF (norm1.LT.d) THEN
1778
                       I1=Jat
1779
                       I0=Iat
1780
                       d=Norm1
1781
                    END IF
1782
                 END IF
1783
              END DO
1784
           END DO
1785

    
1786
           Iat=I0
1787
           n1=I1
1788

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

    
1791

    
1792
           SELECT CASE(idx_zmat(n1))
1793
           CASE (1)
1794
              n2=ind_zmat(2,1)
1795
              n3=ind_zmat(3,1)
1796
           CASE (2)
1797
              n2=ind_zmat(1,1)
1798
              n3=ind_zmat(3,1)
1799
           CASE (3:)
1800
              n2=ind_zmat(idx_zmat(n1),2)
1801
              n3=ind_zmat(idx_zmat(n1),3)
1802
           END SELECT
1803

    
1804
           izm=izm+1
1805
           !        ind_zmat(izm,1)=iat
1806
           !        ind_zmat(izm,2)=n1
1807
           !        ind_zmat(izm,3)=n2
1808
           !        ind_zmat(izm,4)=n3
1809
           Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
1810
           idx_zmat(iat)=izm
1811

    
1812

    
1813
           DejaFait(iat)=.TRUE.
1814
           IdxCaFaire=2
1815
           CaFaire(1)=iat
1816
           CaFaire(2)=0
1817
           IaFaire=1
1818
           FCaf(Iat)=.TRUE.
1819
           DO WHILE (CaFaire(IaFaire).NE.0)
1820
              n1=CaFaire(IaFaire)
1821
              if (debug) WRITE(*,*) 'Iafaire,n1',IaFaire,n1, &
1822
                   DejaFait(n1),idx_zmat(n1)
1823
              if (debug) WRITE(*,*) 'Cafaire', &
1824
                   (CaFaire(J),J=Iafaire,IdxCAfaire)
1825
              SELECT CASE(idx_zmat(n1))
1826
              CASE (1)
1827
                 n2=ind_zmat(2,1)
1828
                 n3=ind_zmat(3,1)
1829
              CASE (2)
1830
                 n2=ind_zmat(1,1)
1831
                 n3=ind_zmat(3,1)
1832
              CASE (3:)
1833
                 n2=ind_zmat(idx_zmat(n1),2)
1834
                 n3=ind_zmat(idx_zmat(n1),3)
1835
              END SELECT
1836

    
1837

    
1838
              DO I3=1,Liaisons(n1,0)
1839
                 IAt=Liaisons(n1,I3)
1840
                 ! PFL 2007.Apr.16 ->
1841
                 ! We add ALL atoms linked to n1 to CaFaire
1842
                 ! But we delete their link to n1
1843
!! PFL 2007.Aug.28 ->
1844
! re-add the test on FCaf in order not to put the same atom more than once in 
1845
! CaFaire array.
1846
                 IF (.NOT.FCaf(Iat)) THEN
1847
                    CaFaire(IdxCaFaire)=Iat
1848
                    IdxCaFaire=IdxCaFaire+1
1849
                    CaFaire(IdxCaFaire)=0
1850
                    FCaf(Iat)=.TRUE.
1851
                 END IF
1852
!! <-- PFL 2007.Aug.28
1853

    
1854
                 Call Annul(Liaisons,Iat,n1)
1855
                 Liaisons(iat,0)=Liaisons(Iat,0)-1
1856
                 if (debug) WRITE(*,*) 'Adding ',iat,' to CaFaire. I3=',I3
1857
                 if (debug) WRITE(*,*) 'Cafaire',(CaFaire(J),J=Iafaire,IdxCAfaire)
1858

    
1859
                 ! <- PFL 2007.Apr.16
1860
                 IF (.NOT.DejaFait(Liaisons(n1,i3))) THEN
1861
                    IAt=Liaisons(n1,I3)
1862
                    izm=izm+1
1863
                    !                 ind_zmat(izm,1)=iat
1864
                    !                 ind_zmat(izm,2)=n1
1865
                    !                 ind_zmat(izm,3)=n2
1866
                    !                 ind_zmat(izm,4)=n3
1867
                    Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
1868
                    idx_zmat(iat)=izm
1869

    
1870
                    n3=iat
1871
                    DejaFait(Iat)=.TRUE.
1872
                 END IF
1873

    
1874
              if (debug) WRITE(*,*) 'I3,n1,Liaisons(n1,:)', &
1875
                   I3,n1,(Liaisons(n1,J),J=0,Liaisons(n1,0))
1876
              END DO
1877
              IaFaire=IaFaire+1
1878

    
1879
              if (debug) THEN
1880
                 WRITE(*,*) 'ind_zmat 7', izm
1881
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
1882
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
1883
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2),ind_zmat(3,3)
1884
                 DO Ip=4,izm
1885
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1), &
1886
                         ind_zmat(Ip,2),  &
1887
                         ind_zmat(Ip,3), ind_zmat(Ip,4)
1888
                 END DO
1889
              END IF
1890

    
1891
           END Do              ! DO WHILE CaFaire
1892
        END IF                 ! FrozBlock(IFrag,0).EQ.0
1893

    
1894
     END DO                    ! Loop on all fragments without frozen atoms
1895

    
1896
     if (debug) THEN
1897
        WRITE(*,*) Na+Izm
1898
        WRITE(*,*) 'Done... ', izm
1899
        DO J=1,Na
1900
           WRITE(*,'(1X,A5,3(1X,F10.6))') Nom(Atome(J)),X(J),Y(J),Z(J)
1901
        END DO
1902
        DO I=1,Izm
1903
           J=ind_zmat(I,1)
1904
           WRITE(*,'(1X,A5,3(1X,F10.6))') 'X ',X(J),Y(J),Z(J)
1905
        END DO
1906
        IF (izm.NE.Na) THEN
1907
           WRITE(*,*) "Atoms not done:"
1908
           DO I=1,Na
1909
              IF (.NOT.DejaFait(I)) WRITE(*,'(I5)',ADVANCE='NO') I
1910
           END DO
1911
        END IF
1912
     END IF
1913

    
1914

    
1915
     ! We have ind_zmat. We calculate val_zmat :-)
1916
     if (debug) WRITE(*,*) "Calculating val_zmat"
1917

    
1918
     val_zmat(1,1)=0.d0
1919
     val_zmat(1,2)=0.d0
1920
     val_zmat(1,3)=0.d0
1921
     val_zmat(2,2)=0.d0
1922
     val_zmat(2,3)=0.d0
1923
     val_zmat(3,3)=0.d0
1924

    
1925
     n1=ind_zmat(2,1)
1926
     n2=ind_zmat(2,2)
1927

    
1928
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
1929

    
1930
     val_zmat(2,1)=norm1
1931

    
1932

    
1933
     n1=ind_zmat(3,1)
1934
     n2=ind_zmat(3,2)
1935
     n3=ind_zmat(3,3)
1936

    
1937
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
1938

    
1939
     val_zmat(3,1)=norm1
1940

    
1941
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
1942
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
1943

    
1944
     val_zmat(3,2)=val
1945

    
1946
     DO i=4,na
1947

    
1948
        n1=ind_zmat(i,1)
1949
        n2=ind_zmat(i,2)
1950
        n3=ind_zmat(i,3)
1951
        n4=ind_zmat(i,4)
1952

    
1953
        if (debug) WRITE(*,*) "Doing i,n1,n2,n3,n4",i,n1,n2,n3,n4
1954
        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
1955

    
1956
        CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
1957
        val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
1958

    
1959
        CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
1960
        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
1961
             vx4,vy4,vz4,norm4)
1962
        CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
1963
             vx5,vy5,vz5,norm5)
1964

    
1965
        val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &
1966
             vx2,vy2,vz2,norm2)
1967

    
1968
        !     write(*,11) n1,n2,norm1,n3,val,n4,val_d
1969
        !11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
1970

    
1971
        val_zmat(i,1)=norm1
1972
        val_zmat(i,2)=val
1973
        val_zmat(i,3)=val_d
1974

    
1975
     END DO
1976

    
1977
     if (debug) THEN
1978
        WRITE(*,*) 'DBG Cre_Zmat_Constr: ind_zmat'
1979
        DO I=1,na
1980
           WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
1981
        END DO
1982

    
1983
        WRITE(*,*) 'DBG Cre_Zmat_Constr: Full zmat'
1984
        DO I=1,na
1985
           WRITE(*,'(1X,I5,1X,I5,F8.4,2(1X,I5,1X,F7.2))') ind_zmat(i,1),(ind_zmat(i,j+1),val_zmat(i,j),j=1,3)
1986
        END DO
1987

    
1988
     END IF
1989

    
1990
     if (debug) WRITE(*,*) "Deallocate (FrozDist,Fragment, NbAtFrag,FragAt)"
1991
     DEALLOCATE(FrozDist,Fragment, NbAtFrag,FragAt)
1992
     if (debug) WRITE(*,*) "Deallocate FrozFrag, IdxFragAt, FrozBlock"
1993
     !  DEALLOCATE(FrozFrag, IdxFragAt, FrozBlock)
1994
     DEALLOCATE(FrozFrag,FrozBlock)
1995
     if (debug) WRITE(*,*) "Deallocate (DistFroz,Liaisons)"
1996
     DEALLOCATE(DistFroz,Liaisons)
1997
     if (debug) WRITE(*,*) "Deallocate (LiaisonsIni)"
1998
     DEALLOCATE(LiaisonsIni)
1999
     if (debug) WRITE(*,*) "Deallocate(CaFaire,DejaFait,FrozAt)"
2000
     DEALLOCATE(CaFaire,DejaFait,FrozAt)
2001
     if (debug) WRITE(*,*) "Deallocate (LiaisonsBis"
2002
     DEALLOCATE(LiaisonsBis)
2003

    
2004

    
2005
     if (debug) WRITE(*,*) "=============================== Exiting Calc_zmat_contr_frag ========================"
2006

    
2007
   END SUBROUTINE Calc_Zmat_const_frag
2008

    
2009