Statistiques
| Révision :

root / src / Calc_zmat_frag.f90 @ 8

Historique | Voir | Annoter | Télécharger (40,61 ko)

1
SUBROUTINE Calc_Zmat_frag(na,atome,ind_zmat,val_zmat,x,y,z,r_cov,fact)
2

    
3
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
4
!
5
!   Goal: to compute a viable Z-Matrix starting from the 
6
!         cartesian coordinates of the atoms
7
!
8
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
9
!
10
! Input:
11
!  na    : Number or atoms 
12
! atome  : Mass number of the atoms (H=1, He=2,...)
13
! x,y,z  : cartesian coordinates of the atoms
14
! r_cov  : array containing the covalent radii of the atoms
15
!  fact  : multiplicative factor used to determine if two atoms are linked.
16
!          see CalcCnct for more details.
17
!
18
! Output:
19
! ind_zmat : INTEGER(na,5) contains the indices of the Z-matrix
20
! val_zmat : REAL(na,3)    contains the values of the Z-matrix
21
!
22
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23
! History
24
!
25
! v1.0 written for Cart a long time ago. Does not use fragment analysis, but was quite good !
26
!
27
! v2.0 12/2007
28
! Comes from Calc_zmat_constr_frag, based on a fragment analysis of the 
29
! system under study.
30
! It should be more flexible and robust.
31
!
32
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33

    
34
  Use Path_module, only : max_Z, NMaxL, Nom, Pi
35
  Use Io_module
36

    
37
  IMPLICIT NONE
38

    
39
  integer(KINT), INTENT(IN) ::  na,atome(na)
40
  real(KREAL), INTENT(IN)   ::  x(Na),y(Na),z(Na),fact
41
  real(KREAL), INTENT(IN)   ::  r_cov(0:Max_Z)
42

    
43
  INTEGER(KINT), INTENT(OUT) :: ind_zmat(Na,5)
44
  real(KREAL), INTENT(OUT)   :: val_zmat(Na,3)
45

    
46
  INTEGER(KINT) :: idx_zmat(NA)
47
! Frozen contains the indices of frozen atoms
48
! INTEGER(KINT) Frozen(*),NFroz
49
! Number of fragment, Index of the current fragment for loops
50
  INTEGER(KINT) :: NbFrag
51
! Fragment(I) contains the fragment index to which I belongs
52
! NbAtFrag(j) contains the number of atoms of fragment j
53
  INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !na
54
! FragAt(i,:) lists the atoms of fragment i
55
  INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !na,na
56
! MaxLFrag(i,1) contains the maximum of links for an atom for fragment i
57
! MaxLFrag(i,2) is the atom that has this number of linked atoms
58
  INTEGER(KINT), ALLOCATABLE :: MaxLFrag(:,:) !na,2
59
! INTEGER(KINT), ALLOCATABLE :: IdxFragAt(:) !na
60
! INTEGER(KINT), ALLOCATABLE :: FrozBlock(:,:) !(na,0:na)
61
! DistFrag contains the distance between a given atom and some other atoms
62
  REAL(KREAL), ALLOCATABLE :: DistFrag(:) !na
63
! FragDist(I) contains the index of the atoms corresponding to DistFrag(I)
64
  INTEGER(KINT), ALLOCATABLE :: FragDist(:) ! na
65

    
66
  INTEGER(KINT) :: IdxCaFaire, IAfaire
67
  INTEGER(KINT), ALLOCATABLE ::  LIAISONS(:,:) ! (Na,0:NMaxL)
68
! INTEGER(KINT), ALLOCATABLE ::  LiaisonsIni(:,:) ! (Na,0:NMaxL)
69
  INTEGER(KINT), ALLOCATABLE :: CAFaire(:) ! (Na)
70

    
71
  real(KREAL) ::  vx1,vy1,vz1,norm1
72
  real(KREAL) ::  vx2,vy2,vz2,norm2
73
  real(KREAL) ::  vx3,vy3,vz3,norm3
74
  real(KREAL) ::  vx4,vy4,vz4,norm4
75
  real(KREAL) ::  vx5,vy5,vz5,norm5
76
  real(KREAL) ::  val,val_d, d12, d13,d23,d
77
  Logical Debug, FirstAt, DebugGaussian
78
  LOGICAL, ALLOCATABLE :: DejaFait(:), FCaf(:) !(na)
79
!  Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen
80
  LOGICAL F1213, F1223,F1323
81

    
82
  INTEGER(KINT) :: I,J,n0,n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax
83
  INTEGER(KINT) :: I3, I1, Ip
84
  INTEGER(KINT) :: I0, Izm, JAt
85
  INTEGER(KINT) :: OrderZmat
86

    
87
  REAL(KREAL) :: sDihe
88

    
89
  INTERFACE
90
     function valid(string) result (isValid)
91
       CHARACTER(*), intent(in) :: string
92
       logical                  :: isValid
93
     END function VALID
94

    
95
     FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
96
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
97
       real(KREAL) ::  v1x,v1y,v1z,norm1
98
       real(KREAL) ::  v2x,v2y,v2z,norm2
99
       real(KREAL) ::  angle
100
     END FUNCTION ANGLE
101

    
102
     FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
103
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
104
       real(KREAL) ::  v1x,v1y,v1z,norm1
105
       real(KREAL) ::  v2x,v2y,v2z,norm2
106
       real(KREAL) ::  SinAngle
107
     END FUNCTION SINANGLE
108

    
109
     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
110
       INTEGER, PARAMETER :: KREAL=KIND(1.0D0)
111
       real(KREAL) ::  v1x,v1y,v1z,norm1
112
       real(KREAL) ::  v2x,v2y,v2z,norm2
113
       real(KREAL) ::  v3x,v3y,v3z,norm3
114
       real(KREAL) ::  angle_d,ca,sa
115
     END FUNCTION ANGLE_D
116
  END INTERFACE
117

    
118
  debug=valid("calc_zmat").OR.valid("calc_zmat_frag")
119
  debugGaussian=valid("zmat_gaussian")
120
  
121
  if (debug) WRITE(*,*) "=============================== Entering Calc_zmat_frag ========================"
122
  if (na.le.2) THEN
123
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
124
     RETURN
125
  END IF
126

    
127
  ALLOCATE(FragDist(Na),Fragment(na), NbAtFrag(na),FragAt(na,na))
128
! ALLOCATE(FrozFrag(na,3), IdxFragAt(na), FrozBlock(na,0:na))
129
! ALLOCATE(FrozFrag(na,3))
130
  ALLOCATE(DistFrag(na),Liaisons(na,0:NMaxL))
131
! ALLOCATE(Liaisons(na,0:NMaxL),LiaisonsIni(na,0:NMaxL))
132
! ALLOCATE(CaFaire(na),DejaFait(Na),FCaf(na),FrozAt(na))
133
  ALLOCATE(CaFaire(na+1),DejaFait(Na),FCaf(na))
134

    
135
  if (debug) THEN
136
     WRITE(*,*) "DBG Calc_zmat_frag - Cartesian coordinates"
137
     DO I=1,na
138
        WRITE(*,'(1X,A3,3(1X,F15.8))') Nom(atome(i)),x(i),y(i),z(i)
139
     END DO
140
  END if
141

    
142
  DO I=1,na
143
     DO J=1,5
144
        Ind_Zmat(I,J)=0
145
     END DO
146
  END DO
147
  
148
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
149
!
150
!     Easy case : 3 or less atoms
151
!
152
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
153
  if (Na.eq.3) THEN
154
     d12=sqrt((x(1)-x(2))**2+(y(1)-y(2))**2+(z(1)-z(2))**2)
155
     d13=sqrt((x(1)-x(3))**2+(y(1)-y(3))**2+(z(1)-z(3))**2)
156
     d23=sqrt((x(3)-x(2))**2+(y(3)-y(2))**2+(z(3)-z(2))**2)
157
     F1213=(d12<=d13)
158
     F1323=(d13<=d23)
159
     F1223=(d12<=d23)
160
     if (debug) THEN
161
        WRITE(*,*) "DEBUG Calc_Zmat 3 atoms"
162
        WRITE(*,*) "d12,d13,d23:",d12,d13,d23
163
        WRITE(*,*) "F1213,F1323,F1223:",F1213,F1323,F1223
164
     END IF
165
     OrderZmat=0
166
     if (F1213) orderZmat=OrderZmat+4
167
     if (F1323) orderZmat=OrderZmat+2
168
     if (F1223) orderZmat=OrderZmat+1
169
     if (debug) WRITE(*,*) 'OrderZmat=',OrderZmat
170
     SELECT CASE (OrderZmat)
171
        CASE (0)
172
! F F F ordre 2-3----1 
173
           ind_zmat(1,1)=3
174
           ind_zmat(2,1)=2
175
           ind_zmat(2,2)=3
176
           ind_zmat(3,1)=1
177
           ind_zmat(3,2)=3
178
           ind_zmat(3,3)=2
179
        CASE (2)
180
! F T F ordre 1-3----2
181
           ind_zmat(1,1)=3
182
           ind_zmat(2,1)=1
183
           ind_zmat(2,2)=3
184
           ind_zmat(3,1)=2
185
           ind_zmat(3,2)=3
186
           ind_zmat(3,3)=1
187
        CASE (3)
188
! F T T ordre 2---1-3
189
           ind_zmat(1,1)=1
190
           ind_zmat(2,1)=3
191
           ind_zmat(2,2)=1
192
           ind_zmat(3,1)=2
193
           ind_zmat(3,2)=1
194
           ind_zmat(3,3)=3
195
        CASE (5)
196
! T F T  ordre 1-2----3
197
           ind_zmat(1,1)=2
198
           ind_zmat(2,1)=1
199
           ind_zmat(2,2)=2
200
           ind_zmat(3,1)=3
201
           ind_zmat(3,2)=2
202
           ind_zmat(3,3)=1
203
        CASE (7)
204
! T T T ordre 3----1-2
205
           ind_zmat(1,1)=1
206
           ind_zmat(2,1)=2
207
           ind_zmat(2,2)=1
208
           ind_zmat(3,1)=3
209
           ind_zmat(3,2)=1
210
           ind_zmat(3,3)=2
211
        END SELECT
212

    
213
     IF (debug) THEN
214
        WRITE(*,*) "DBG Calc_zmat_frag - Nat=3 -"
215
        DO i=1,na
216
           WRITE(*,'(1X,A3,5(1X,I5))') Nom(Atome(ind_zmat(i,1))),(ind_zmat(i,j),j=1,5)
217
        END DO
218
     END IF
219

    
220
     ! We have ind_zmat, we fill val_zmat
221
     val_zmat=0.d0
222
     n2=ind_zmat(2,1)
223
     n1=ind_zmat(2,2)
224
     d=sqrt((x(n1)-x(n2))**2+(y(n1)-y(n2))**2+(z(n1)-z(n2))**2)
225
     val_zmat(2,1)=d
226
     n1=ind_zmat(3,1)
227
     n2=ind_zmat(3,2)
228
     n3=ind_zmat(3,3)
229
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
230
     if (debug) write(*,*) n1,n2,norm1
231
     val_zmat(3,1)=norm1
232

    
233
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
234
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
235
     if (debug) write(*,*) n2,n3,norm2,val
236
     val_zmat(3,2)=val
237

    
238
     RETURN
239
  END IF !matches if (Na.eq.3) THEN
240
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
241
!
242
!  End of Easy case : 3 or less atoms
243
!
244
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
245
!
246
! General Case 
247
!
248
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
249
!
250

    
251
! Initialization
252
  DejaFait=.False.
253
  Liaisons=0
254
  ind_zmat=0
255
  val_zmat=0.d0
256

    
257
  if (debug) WRITE(*,*) "Coucou from Calc_zmat_frag.f90; L240"
258
  
259
  if (debug) THEN
260
     WRITE(*,*) "Bonds initialized"
261
     WRITE(*,*) 'Covalent radii used'
262
     DO iat=1,na
263
        i=atome(iat)
264
        WRITE(*,*) Nom(I),Iat,r_cov(i),r_cov(i)*fact
265
     END DO
266
  END IF
267

    
268
1003 FORMAT(1X,I4,' - ',25(I5))
269

    
270
!   First step : connectivity  are calculated
271

    
272
  Call CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)
273

    
274
  if (debug) THEN
275
     WRITE(*,*) "Connections calculated"
276
     DO IL=1,na
277
        WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
278
     END DO
279
  END IF
280

    
281
  FCaf=.TRUE.
282
  Call Decomp_frag(na,liaisons,FCaf,nbfrag,Fragment,NbAtFrag,FragAt)
283
  
284
  IF (debug) THEN
285
     WRITE(*,*) 'I found ',NbFrag, 'fragments'
286
     WRITE(*,*) (NbAtFrag(I),I=1,NbFrag)
287
     DO I=1,NbFrag
288
        WRITE(*,*) NbAtFrag(I)
289
        WRITE(*,*) 'Fragment ', i
290
        DO J=1,Na
291
           IF (Fragment(J).EQ.I) WRITE(*,'(1X,I3,1X,A5,3(1X,F10.6))') J,Nom(Atome(J)) &
292
                ,X(J),Y(J),Z(J)
293
        END DO
294
        WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
295
     END DO
296
  END IF
297

    
298
 ALLOCATE(MaxLFrag(NbFrag,2))
299

    
300
 MaxLFrag=0
301

    
302
  DO I=1,NbFrag
303
     MaxLFrag(I,1)=Liaisons(FragAt(I,1),0)
304
     MaxLFrag(I,2)=FragAt(I,1)
305

    
306
     DO J=1,NbAtFrag(I)
307
      Iat=FragAt(I,J)
308
        IF (Liaisons(IAt,0).GT.MaxLFrag(I,1)) THEN
309
       MaxLFrag(I,1)=Liaisons(Iat,0)
310
       MaxLFrag(I,2)=Iat
311
    END IF
312
     END DO
313
     IF (debug) WRITE(*,*) 'Frag :',I,', atom ',MaxLFrag(I,2), ' has ',MaxLFrag(I,2),' links'
314
  END DO
315

    
316
  !     We will now build the molecule fragment by fragment
317
  !     We choose the starting fragment with two criteria:
318
  !     1- Number of linked atoms:
319
  !       * >=3 is good as it fully defines the coordinate space
320
  !       * 2 is ok as we can either use a 3rd atom from the same fragment
321
  !       or add a X atom somewhere but this complicates quite a lot the way
322
  !       to treat the conversion from cartesian to zmat latter
323
  !       * 1 is bad...
324
  !     2- Size of the fragment
325
  !     this allows us to deal more easily with cases 1- when number of 
326
  !     directly linked atoms is less than 3
327

    
328
  IFrag=0
329
! I0 is the number of connections of the best fragment
330
  I0=0
331
! I1 is the number of atoms of the best fragment
332
  I1=0
333
  IAt=0
334
  DO I=1,NbFrag
335
    SELECT CASE(MaxLFrag(I,1)-I0)
336
     CASE (1:)
337
        IFrag=I
338
        I0=MaxLFrag(I,1)
339
        I1=NbAtFrag(I)
340
        IAt=MaxLFrag(I,2)
341
     CASE (0)
342
        if (NbAtFrag(I).GT.I1) THEN
343
           IFrag=I
344
           I0=MaxLFrag(I,1)
345
           I1=NbAtFrag(I)
346
           IAt=MaxLFrag(I,2)
347
        END IF
348
     END SELECT
349
  END DO
350

    
351
  if (debug) WRITE(*,'(1X,A,I5,A,I5,A,I5,A,I5)') 'Starting with fragment:',IFrag,' with ',I0   &
352
      ,' max links for atom',IAt,' fragment size',NbAtFrag(IFrag)
353

    
354
  !     We will build the first fragment in a special way, as it will
355
  !     set the coordinates system
356

    
357
  if (debug) WRITE(*,*) 'Fragment 1, starting with atom:',IAt, &   
358
    'with ',I0,' connections'
359

    
360
  DejaFait=.FALSE.
361
  FCaf=.FALSE.
362
  
363
  izm=0
364
  SELECT CASE (I0)
365
  CASE(3:)
366
     if (debug) WRITE(*,*) 'DBG select case I0 3'
367
     n0=Iat
368

    
369
     ITmp=2
370
     sDihe=0.
371
     n2=IAt
372
     n3=Liaisons(Iat,1)
373
     ! We search for the third atom while making sure that it is not aligned with first two !
374
     DO While ((ITmp.LE.Liaisons(Iat,0)).AND.(sDihe.LE.0.09d0))
375
        n4=Liaisons(Iat,Itmp)
376
        CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
377
        CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
378
        val_d=angle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
379
        sDiHe=abs(sin(val_d*pi/180.d0))
380
        if (debug) Write(*,*) 'Trying n2,n3,n4,sdihe,val_d',n2,n3,n4,sdihe,val_d
381
        Itmp=Itmp+1
382
     END DO
383
     If (debug) WRITE(*,*) 'Itmp,Liaisons',Itmp,Liaisons(Iat,1:NMaxL)
384
     Liaisons(Iat,Itmp-1)=Liaisons(iat,2)
385
     Liaisons(Iat,2)=n4
386
     If (debug) WRITE(*,*) 'Itmp,Liaisons',Itmp,Liaisons(Iat,1:NMaxL)
387

    
388
     if (sDihe.LE.0.09d0) THEN
389
        WRITE(*,*) "PROBLEM !!! All atoms linked to ",n0," are aligned..."
390
    WRITE(*,*) "Surprising as this atom has at least three bonds... For NOW STOP"
391
    STOP
392
     END IF
393

    
394
     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,  vx5,vy5,vz5,norm5)
395

    
396

    
397
! We search for the fourth atom while checking that it is not aligned with 1st and 2nd atoms.
398
     Itmp=3
399
     sDiHe=0.
400
! PFL 28 Dec 2007 ->
401
! I had a test on the dihedral angle, but I cannot see why it is important to have
402
! a non planar fragment at the begining... ethylene works and is fully planar
403
! I thus suppress this test
404
!
405
!     DO While ((ITmp.LE.Liaisons(Iat,0)).AND.(sDihe.LE.0.09d0))
406
!        ITmp=ITmp+1
407
        n1=Liaisons(Iat,Itmp)
408
        if (debug) WRITe(*,*) 'trying n1,n2,n3,n4',n1,n2,n3,n4
409
        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
410
        ! Is this atom aligned with n2-n3 ?
411
        val_d=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
412
        sDiHe=abs(sin(val_d*pi/180.d0))
413
        if (debug) Write(*,*) 'Angle n3-n2-n1',val_d
414
        if (sDiHe.le.0.09d0) THEN
415
           ! As atoms n2,n3 and n4 are not aligned, we interchange n3 and n4 so that n1,n2 and n3 are not aligned
416
           if (debug) WRITE(*,*) "n3-n2-n1 aligned, we interchange n3 and n4"
417
           n1=n3
418
           n3=n4
419
           n4=n1
420
           n1=Liaisons(Iat,ITmp)
421
           CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)           
422
           CALL vecteur(n2,n4,x,y,z,vx3,vy3,vz3,norm3)           
423
           val_d=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
424
           if (debug) Write(*,*) 'NEW Angle n3-n2-n1',val_d
425
        END IF
426

    
427
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
428
        ! avant de tester l'angle diedre, il faut verifier que ce 4e atome n'est pas
429
        ! aligne avec les 2 premiers. 
430
        ! comme on a deja test? que les 3 premiers ne sont pas alignes,
431
        ! s'il est align? avec les 2 premiers, on peut inverser le role de 2 et 3.
432
        ! On pourrait tout simplier en faisant une bete recherche parmi tous les atomes geles
433
        ! de ce bloc (au moins 4 ?) avec le critere 1) on les range par distance croissante
434
        ! 2) on les scanne tant que l'angle valence n'est pas bon, puis tant que diehedre pas bon
435
        ! on pourrait comme ca faire un tableau avec les atomes ranges d'abord pour le fragment s?lectionn?
436
        ! puis les atomes des autres fragment par distance croissante.
437
        ! les autres fragments ne seraient additonn?s que si l'on ne trouve pas notre bonheur dans le premier bloc
438
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
439

    
440
        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
441
        val_d=angle_d(vx4,vy4,vz4,norm4,vx5,vy5,vz5,norm5, &
442
             vx2,vy2,vz2,norm2)
443
        sDihe=abs(sin(val_d*pi/180.d0))
444
        if (debug) WRITE(*,*) 'n2,n3,n4,n1, angle_d',n2,n3,n4,n1,val_d
445
!     END DO
446

    
447
     DejaFait(n2)=.TRUE.
448
     DejaFait(n3)=.TRUE.
449
     DejaFait(n4)=.TRUE.
450

    
451
!     if (sDihe.LE.0.09d0) THEN
452
!        !     None of the atoms linked to IAt are good to define the third
453
!        !     direction in space...
454
!        !     We will look at the other atoms
455
!        !     we might improve the search so as to take the atom closest to IAt
456
!        if (debug) WRITE(*,*) "All atoms linked to ",Iat," are in a plane. Looking for other atoms"
457
!        ITmp=0
458
!        DO I=1,NbAtFrag(IFrag)
459
!           JAt=FragAt(IFrag,I)
460
!           if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
461
!              n1=JAt
462
!              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
463
!              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
464
!              val_d=angle_d(vx4,vy4,vz4,norm4, &
465
!                   vx5,vy5,vz5,norm5, &
466
!                   vx2,vy2,vz2,norm2)
467
!              if (abs(sin(val_d)).GE.0.09D0) THEN
468
!                 ITmp=ITmp+1
469
!                 DistFrag(ITmp)=Norm1
470
!                 FragDist(ITmp)=JAt
471
!              END IF
472
!           END IF
473
!        END DO
474

    
475
!        IF (ITmp.EQ.0) THEN
476
!           !     All dihedral angles between frozen atoms are less than 5?
477
!           !     (or more than 175?). We have to look at other fragments :-(
478
!           DO I=1,NFroz
479
!              Jat=Frozen(I)
480
!              if (.NOT.DejaFait(JAt)) THEN
481
!                 n1=JAt
482
!                 CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
483
!                 CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
484
!                 val_d=angle_d(vx4,vy4,vz4,norm4, &
485
!                      vx5,vy5,vz5,norm5, &
486
!                      vx2,vy2,vz2,norm2)
487
!                 if (abs(sin(val_d)).GE.0.09D0) THEN
488
!                    ITmp=ITmp+1
489
!                    DistFrag(ITmp)=Norm1
490
!                    FragDist(ITmp)=JAt
491
!                 END IF
492
!              END IF
493
!           END DO
494
!           IF (ITmp.EQ.0) THEN
495
!              !     All frozen atoms are in a plane... too bad
496
!              WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
497
!              WRITE(*,*) 'For now, I do not treat this case'
498
!              STOP
499
!           END IF
500
!        END IF
501
!        I1=0
502
!        d=1e9
503
!        DO I=1,ITmp
504
!           IF (DistFrag(I).LE.d) THEN
505
!              d=DistFrag(I)
506
!              I1=FragDist(I)
507
!           END IF
508
!        END DO
509
!     ELSE                !(sDihe.LE.0.09d0)
510
!        I1=FrozBlock(IAt,ITmp)
511
!        if (debug) WRITE(*,*) 'I1,n1:',I1,n1
512
!     END IF                 !(sDihe.LE.0.09d0)
513
!     !     we now have the atom that is closer to the first one and that
514
!     !     forms a non 0/Pi dihedral angle
515
!
516
!  <- PFL 28 Dec 2007
517

    
518
! We construct the begining of the Z-Matrix
519

    
520
     ind_zmat(1,1)=n2
521
     ind_zmat(2,1)=n3
522
     ind_zmat(2,2)=n2
523
     ind_zmat(3,1)=n4
524
     ind_zmat(3,2)=n2
525
     ind_zmat(3,3)=n3
526
     DejaFait(n2)=.TRUE.
527
     DejaFait(n3)=.TRUE.
528
     DejaFait(n4)=.TRUE.
529
     CaFaire(1)=n3
530
     CaFaire(2)=n4
531

    
532
! PFL 28 Dec 2007
533
! We have not selected a fourth atom, so that the following is not needed
534
!     ind_zmat(4,1)=I1
535
!     ind_zmat(4,2)=n2
536
!     ind_zmat(4,3)=n3
537
!     ind_zmat(4,4)=n4
538
!     DejaFait(I1)=.TRUE.
539
!     CaFaire(3)=I1
540
!     CaFaire(4)=0
541
!     IdxCaFaire=4
542
!     izm=4
543
!     FCaf(I1)=.TRUE.
544
!!!!!!!
545
! and replaced by:
546
     CaFaire(3)=0
547
   IdxCaFaire=3
548
   izm=3
549
!
550
! <- PFL 28 Dec 2007
551
   
552
     FCaf(n2)=.TRUE.
553
     FCaf(n3)=.TRUE.
554
   FirstAt=.TRUE.
555
     DO I=3,Liaisons(Iat,0)
556
        IF (.NOT.DejaFait(Liaisons(Iat,I))) THEN
557
           izm=izm+1
558
           !           ind_zmat(izm,1)=Liaisons(Iat,I)
559
           !           ind_zmat(izm,2)=n2
560
           !           ind_zmat(izm,3)=n3
561
           !           ind_zmat(izm,4)=n4
562
           Call add2indzmat(na,izm,Liaisons(Iat,I),n2,n3,n4,ind_zmat,x,y,z)
563
           if (FirstAt) THEN
564
           n4=Liaisons(Iat,I)
565
         FirstAt=.FALSE.
566
       END IF
567
           IF (.NOT.FCaf(Liaisons(Iat,I))) THEN
568
              CaFaire(IdxCaFaire)=Liaisons(Iat,I)
569
              IdxCaFaire=IdxCaFaire+1
570
              CaFaire(IdxCaFaire)=0
571
              FCaf(Liaisons(Iat,I))=.TRUE.
572
           END IF
573
           DejaFait(Liaisons(Iat,I))=.TRUE.
574
        END IF
575
     END DO
576

    
577
     if (debug) THEN
578
        WRITE(*,*) "Ind_zmat 0 - SELECT CASE I0 3: -- izm=",izm
579
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
580
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
581
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
582
        DO I=4,izm
583
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2),  &
584
                ind_zmat(I,3), ind_zmat(I,4)
585
        END DO
586
     END IF
587

    
588

    
589
     !     First four atoms (at least) have been put we go on with next parts
590
     !     of this fragment... later
591

    
592

    
593
  CASE(2)
594
     if (debug) WRITE(*,*) 'DBG select case I0 2'
595
   WRITE(*,*) "PFL 28 Dec 2007: No test of alignment here :(  TO DO TO DO TO DO"
596
     ind_zmat(1,1)=IAt
597
     ind_zmat(2,1)=Liaisons(IAt,1)
598
     ind_zmat(2,2)=IAt
599
     ind_zmat(3,1)=Liaisons(IAt,2)
600
     ind_zmat(3,2)=IAt
601
     ind_zmat(3,3)=Liaisons(IAt,1)
602
     DejaFait(IAt)=.TRUE.
603
     DejaFait(Liaisons(Iat,1))=.TRUE.
604
     DejaFait(Liaisons(Iat,2))=.TRUE.
605
     CaFaire(1)=Liaisons(Iat,1)
606
     CaFaire(2)=Liaisons(Iat,2)
607
     FCaf(Liaisons(Iat,1))=.TRUE.
608
     FCaf(Liaisons(Iat,2))=.TRUE.
609

    
610
! PFL 28 Dec 2007 ->
611
! We do NOT need a fourth atom !!! The third direction in space is defined by the cross
612
! product of the first two directions
613
!
614
! the following is thus commented
615
!
616
!     !     We search for a fourth atom, first in the FrozBlock atoms
617
!     ITmp=2
618
!     sDihe=0.
619
!     n2=IAt
620
!     n3=Liaisons(Iat,1)
621
!     n4=Liaisons(Iat,2)
622
!     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
623
!     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
624
!     CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2, vx5,vy5,vz5,norm5)
625
!
626
!     DO While ((ITmp.LE.FrozBlock(Iat,0)).AND.(sDihe.LE.0.09d0))
627
!        ITmp=ITmp+1
628
!        n1=FrozBlock(Iat,Itmp)
629
!        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
630
!        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
631
!        val_d=angle_d(vx4,vy4,vz4,norm4,  &
632
!             vx5,vy5,vz5,norm5,           &
633
!             vx2,vy2,vz2,norm2)
634
!        sDihe=abs(sin(val_d))
635
!     END DO
636
!     IF (debug) WRITE(*,*) 'DBG 4th atom, ITmp, sDihe',ITmp, sDihe
637
!     if (sDihe.LE.0.09d0) THEN
638
!        !     None of the frozen atoms linked to IAt are good to define the third
639
!        !     direction in space...
640
!        !     We will look at the other frozen atoms
641
!        !     we might improve the search so as to take the atom closest to IAt
642
!        ITmp=0
643
!        DO I=1,NbAtFrag(IFrag)
644
!           JAt=FragAt(IFrag,I)
645
!           if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
646
!              n1=JAt
647
!              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
648
!              CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
649
!              val_d=angle_d(vx4,vy4,vz4,norm4,    &
650
!                   vx5,vy5,vz5,norm5,              &
651
!                   vx2,vy2,vz2,norm2)
652
!              if (abs(sin(val_d)).GE.0.09D0) THEN
653
!                 ITmp=ITmp+1
654
!                 DistFrag(ITmp)=Norm1
655
!                 FragDist(ITmp)=JAt
656
!              END IF
657
!           END IF
658
!        END DO
659
!        IF (ITmp.EQ.0) THEN
660
!           !     All dihedral angles between frozen atoms are less than 5?
661
!           !     (or more than 175?). We have to look at other fragments :-(
662
!           DO I=1,NFroz
663
!              Jat=Frozen(I)
664
!              if (.NOT.DejaFait(JAt)) THEN
665
!                 n1=JAt
666
!                 CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
667
!                 CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
668
!                 val_d=angle_d(vx4,vy4,vz4,norm4,       &
669
!                      vx5,vy5,vz5,norm5,                 &
670
!                      vx2,vy2,vz2,norm2)
671
!                 if (abs(sin(val_d)).GE.0.09D0) THEN
672
!                    ITmp=ITmp+1
673
!                    DistFrag(ITmp)=Norm1
674
!                    FragDist(ITmp)=JAt
675
!                 END IF
676
!              END IF
677
!           END DO
678
!           IF (ITmp.EQ.0) THEN
679
!              !     All frozen atoms are in a plane... too bad
680
!              WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
681
!              WRITE(*,*) 'For now, I do not treat this case'
682
!              STOP
683
!           END IF
684
!        END IF
685
!        I1=0
686
!        d=1e9
687
!        DO I=1,ITmp
688
!           IF (DistFrag(I).LE.d) THEN
689
!              d=DistFrag(I)
690
!              I1=FragDist(I)
691
!           END IF
692
!        END DO
693
!     ELSE                   !(sDihe.LE.0.09d0)
694
!        I1=FrozBlock(IAt,ITmp)
695
!     END IF                 !(sDihe.LE.0.09d0)
696
!     !     we now have the atom that is closer to the first one and that
697
!     !     forms a non 0/Pi dihedral angle
698
!     !     ind_zmat(4,1)=I1
699
!     !     ind_zmat(4,2)=IAt
700
!     !     ind_zmat(4,3)=Liaisons(Iat,1)
701
!     !     ind_zmat(4,4)=Liaisons(Iat,2)
702
!     n3=Liaisons(Iat,1)
703
!     n4=Liaisons(Iat,2)
704
!     Call add2indzmat(na,4,I1,Iat,n3,n4,ind_zmat,x,y,z)
705
!     Liaisons(Iat,1)=n3
706
!     Liaisons(Iat,2)=n4
707
!     DejaFait(I1)=.TRUE.
708
!     CaFaire(3)=I1
709
!     CaFaire(4)=0
710
!     IdxCaFaire=4
711
!     izm=4
712
!     FCaf(I1)=.TRUE.
713
!
714
!!!!!! <- PFL 28 Dec 2007
715

    
716
     CaFaire(3)=0
717
     IdxCaFaire=3
718
     izm=3
719

    
720
  CASE(1)
721
     if (debug) WRITE(*,*) 'DBG select case I0 1, NbAtFrag=',NbAtFrag(IFrag)
722
     ind_zmat(1,1)=IAt
723
     ind_zmat(2,1)=Liaisons(IAt,1)
724
     ind_zmat(2,2)=IAt
725
     DejaFait(IAt)=.TRUE.
726
     DejaFait(Liaisons(Iat,1))=.TRUE.
727
     IdxCaFaire=2
728
     CaFaire(1)=Liaisons(Iat,1)
729
     CaFaire(2)=0
730
     FCaf(Liaisons(Iat,1))=.TRUE.
731

    
732
! PFL 28 Dec 2007 ->
733
! We do NOT need a fourth atom. So we will look only for a third atom
734
!
735
!!!!
736
!
737
     !     We search for a third and fourth atoms, first in the FrozBlock atoms
738
     !     It should not be possible to have   (FrozBlock(Iat,0).GT.2) and
739
     !     iat linked to only one atom !       
740

    
741

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

    
746
     ITmp=0
747
     CALL vecteur(Liaisons(Iat,1),IAt,x,y,z,vx2,vy2,vz2,norm2)
748

    
749
! PFL 28 Dec 2007: As MaxL=1 I think that there is at most 2 atoms in this fragment... 
750
! so that the following loop is useless... this should be tested more carefully
751
     DO I=1,NbAtFrag(IFrag)
752
        JAt=FragAt(IFrag,I)
753
        if (.NOT.DejaFait(JAt)) THEN
754
           CALL vecteur(JAt,IAt,x,y,z,vx1,vy1,vz1,norm1)
755
           if (abs(cos(angle(vx1,vy1,vz1,norm1,  &
756
                vx2,vy2,vz2,norm2))).LE.0.996d0) THEN
757
              ITmp=ITmp+1
758
              DistFrag(ITmp)=Norm1
759
              FragDist(ITmp)=JAt
760
           END IF
761
        END IF
762
     END DO
763

    
764
     IF (ITMP.EQ.0) THEN
765
        !     We have no atoms correct in this fragment, we use
766
        !     atoms from other fragments
767
        DO Jat=1,Na
768
! DejaFait(Iat)=.TRUE. so that we do not need to test Jat/=Iat
769
           if (.NOT.DejaFait(JAt)) THEN
770
              CALL vecteur(JAt,IAt,x,y,z,vx1,vy1,vz1,norm1)
771
              if (abs(cos(angle(vx1,vy1,vz1,norm1, &
772
                   vx2,vy2,vz2,norm2))).LE.0.996d0) THEN
773
                 ITmp=ITmp+1
774
                 DistFrag(ITmp)=Norm1
775
                 FragDist(ITmp)=JAt
776
              END IF
777
           END IF
778
        END DO
779
        IF (ITMP.EQ.0) THEN
780
           WRITE(*,*) 'It seems all atoms are aligned'
781
           WRITE(*,*) 'Case non treated for now :-( '
782
           STOP
783
        END IF
784
     END IF
785

    
786
     I1=0
787
     d=1e9
788
! PFL 28 Dec 2007: There exists some F90 intrinsics to find the smallest element of an array.
789
! The following loop should be replaced by it !
790
     DO I=1,ITmp
791
        IF (DistFrag(I).LE.d) THEN
792
           I1=FragDist(I)
793
           d=DistFrag(I)
794
        END IF
795
     END DO
796

    
797
     !     we now have the atom that is closer to the first one and that
798
     !     forms a non 0/Pi valence angle
799
     ind_zmat(3,1)=I1
800
     ind_zmat(3,2)=IAt
801
     ind_zmat(3,3)=Liaisons(Iat,1)
802
     DejaFait(I1)=.TRUE.
803
     CaFaire(2)=I1
804
     FCaf(I1)=.TRUE.
805

    
806

    
807
! PFL 28 Dec 2007 ->
808
! We do NOT need a fourth atom so that the following is suppressed
809
!
810
!     !     we search for a fourth atom !
811
!     !     We search for a fourth atom, first in the FrozBlock atoms
812
!     ITmp=2
813
!     sDihe=0.
814
!     n2=IAt
815
!     n3=Liaisons(Iat,1)
816
!     n4=I1
817
!     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
818
!     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
819
!     CALL produit_vect(vx3,vy3,vz3, vx2,vy2,vz2,vx5,vy5,vz5,norm5)
820
!
821
!     !     We will look at the other frozen atoms
822
!     !     we might improve the search so as to take the atom closest to IAt
823
!     ITmp=0
824
!     DO I=1,NbAtFrag(IFrag)
825
!        JAt=FragAt(IFrag,I)
826
!        if (FrozAt(Jat).AND.(.NOT.DejaFait(JAt))) THEN
827
!           n1=JAt
828
!           CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
829
!           CALL produit_vect(vx1,vy1,vz1, vx2,vy2,vz2, vx4,vy4,vz4,norm4)
830
!           val_d=angle_d(vx4,vy4,vz4,norm4,  &
831
!                vx5,vy5,vz5,norm5,   &
832
!                vx2,vy2,vz2,norm2)
833
!           if (abs(sin(val_d)).GE.0.09D0) THEN
834
!              ITmp=ITmp+1
835
!              DistFrag(ITmp)=Norm1
836
!              FragDist(ITmp)=JAt
837
!           END IF
838
!        END IF
839
!     END DO
840
!     IF (ITmp.EQ.0) THEN
841
!        !     All dihedral angles between frozen atoms are less than 5?
842
!        !     (or more than 175?). We have to look at other fragments :-(
843
!        DO I=1,NFroz
844
!           Jat=Frozen(I)
845
!           if (.NOT.DejaFait(JAt)) THEN
846
!              n1=JAt
847
!              CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
848
!              CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, vx4,vy4,vz4,norm4)
849
!              val_d=angle_d(vx4,vy4,vz4,norm4,   &
850
!                   vx5,vy5,vz5,norm5,           &
851
!                   vx2,vy2,vz2,norm2)
852
!              if (abs(sin(val_d)).GE.0.09D0) THEN
853
!                 ITmp=ITmp+1
854
!                 DistFrag(ITmp)=Norm1
855
!                 FragDist(ITmp)=JAt
856
!              END IF
857
!           END IF
858
!       END DO
859
!        IF (ITmp.EQ.0) THEN
860
!           !     All frozen atoms are in a plane... too bad
861
!           WRITE(*,*) 'ERROR: It seems that all frozen atoms are in a plane'
862
!           WRITE(*,*) 'For now, I do not treat this case'
863
!           STOP
864
!        END IF
865
!     END IF                 ! ITmp.EQ.0 after scaning fragment
866
!     I1=0
867
!     d=1e9
868
!     DO I=1,ITmp
869
!        IF (DistFrag(I).LE.d) THEN
870
!           d=DistFrag(I)
871
!           I1=FragDist(I)
872
!        END IF
873
!     END DO
874
!
875
!     !     we now have the atom that is closer to the first one and that
876
!     !     forms a non 0/Pi dihedral angle
877
!     !     ind_zmat(4,1)=I1
878
!     !     ind_zmat(4,2)=IAt
879
!     !     ind_zmat(4,3)=ind_zmat(2,1)
880
!     !     ind_zmat(4,4)=ind_zmat(3,1)
881
!     n3=ind_zmat(2,1)
882
!     n4=ind_zmat(3,1)
883
!     Call add2indzmat(na,4,I1,IAt,n3,n4,ind_zmat,x,y,z)
884
!     ind_zmat(2,1)=n3
885
!     ind_zmat(3,1)=n4
886
!     DejaFait(I1)=.TRUE.
887
!     CaFaire(3)=I1
888
!     CaFaire(4)=0
889
!     IdxCaFaire=4
890
!     izm=4
891
!     FCaf(I1)=.TRUE.
892
!!!!!!!!!!!
893
!
894
! <- PFL 28 Dec 2007
895

    
896
     CaFaire(3)=0
897
   IdxCaFaire=3
898
   
899
  CASE(0)
900
     WRITE(*,*) 'All atoms are separated .. '
901
     WRITE(*,*) 'this case should be treated separately !'
902
     STOP
903
  END SELECT
904

    
905
  if (debug) THEN
906
     WRITE(*,*) 'ind_zmat 1 izm=',izm
907
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
908
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
909
     WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
910
     DO I=4,izm
911
        WRITE(*,'(1X,4(1X,I5))') ind_zmat(I,1), ind_zmat(I,2),  &
912
             ind_zmat(I,3), ind_zmat(I,4)
913
     END DO
914
  END IF
915

    
916
  DO I=1,izm
917
     Idx_zmat(ind_zmat(I,1))=i
918
  END Do
919

    
920
  !     at least first three atoms of this fragment done...
921
  !     we empty the 'cafaire' array before going on
922
  IAFaire=1
923
  DO WHILE (CaFaire(IaFaire).NE.0)
924
     n1=CaFaire(IaFaire)
925
     n2=ind_zmat(idx_zmat(N1),2)
926
     if (idx_zmat(N1).EQ.2) THEN
927
        !     We have a (small) problem: we have to add atoms linked to the 2nd
928
        !     atom of the zmat. This is a pb because we do not know
929
        !     which atom to use to define the dihedral angle
930
        !     we take the third atom of the zmat
931
        n3=ind_zmat(3,1)
932
     ELSE
933
        n3=ind_zmat(idx_zmat(n1),3)
934
     END IF
935
   
936
   FirstAt=.TRUE.
937
     DO I=1,Liaisons(n1,0)
938
        IAt=Liaisons(n1,I)
939
! PFL 29.Aug.2008 ->
940
! We dissociate the test on 'DejaFait' that indicates that this atom
941
! has already been put in the Zmat
942
! from the test on FCaf that indicates that this atom has been put in the
943
! 'CAFaire' list that deals with identifying its connectivity.
944
! Those two test might differ in some cases.
945
        IF (.NOT.DejaFait(Iat)) THEN
946
           izm=izm+1
947
           if (debug) WRITE(*,*) ">1< Adding atom ",Iat,"position izm=",izm
948
           !           ind_zmat(izm,1)=iat
949
           !           ind_zmat(izm,2)=n1
950
           !           ind_zmat(izm,3)=n2
951
           !           ind_zmat(izm,4)=n3
952
           Call add2indzmat(na,izm,iat,n1,n2,n3,ind_zmat,x,y,z)
953
       if (FirstAt) THEN
954
          n3=Iat
955
        FirstAt=.FALSE.
956
           END IF
957
           idx_zmat(iat)=izm
958
           DejaFait(iat)=.TRUE.
959
        END IF
960
        IF (.NOT.FCaf(Iat)) THEN
961
           CaFaire(IdxCaFaire)=iat
962
           IdxCaFaire=IdxCaFaire+1
963
           CaFaire(IdxCaFaire)=0
964
           FCaf(Iat)=.TRUE.
965
        END IF
966
! <- PFL 29.Aug.2008
967
     END DO
968
     IaFaire=IaFaire+1
969
  END Do                    ! DO WHILE CaFaire
970

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

    
982
  !     We have finished putting atoms linked to the first one
983
  ! There should not be any atom left from this fragment. We check:
984
  !     we will add other atoms of this fragment
985
  DO I=1,NbAtFrag(IFrag)
986
     Iat=FragAt(IFrag,I)
987
     if (debug) WRITE(*,*) "DBG: I,Iat,dejafait",I,Iat,DejaFait(Iat)
988
     IF (.NOT.DejaFait(Iat)) THEN
989
           WRITE(*,*) 'Treating atom I,Iat',I,Iat
990
        
991
     END IF
992

    
993
  END DO
994

    
995
  NbAtFrag(Ifrag)=0
996
  MaxLFrag(IFrag,1)=0
997
  
998
  !     we start again with the rest of the molecule...
999
  ! v 1.01 We add the fragment in the order corresponding to NbAtFrag
1000
  KMax=NbFrag-1
1001

    
1002
  IF (DEBUG) WRITE(*,*) "Adding the ",Kmax," remaining fragments"
1003
  DO K=1, KMax
1004
  IFrag=0
1005
  I0=0
1006
  IAt=0
1007
  I1=0
1008
  DO I=1,NbFrag
1009
    SELECT CASE(MaxLFrag(I,1)-I0)
1010
     CASE (1:)
1011
        IFrag=I
1012
        I0=MaxLFrag(I,1)
1013
        IAt=MaxLFrag(I,2)
1014
        I1=NbAtFrag(I)
1015
     CASE (0)
1016
        if (NbAtFrag(I).GT.I1) THEN
1017
           IFrag=I
1018
           I0=MaxLFrag(I,1)
1019
           IAt=MaxLFrag(I,2)
1020
           I1=NbAtFrag(I)
1021
        END IF
1022
     END SELECT
1023
 
1024
  END DO
1025

    
1026
  if (debug) WRITE(*,'(1X,A,I5,A,I5,A,I5,A,I5)') 'Adding fragment:',IFrag,' with ',I0   &
1027
      ,' max links for atom',IAt,' fragment size',NbAtFrag(IFrag)
1028
    
1029
     MaxLFrag(IFrag,1)=0
1030

    
1031
! We search for the closest atoms of the previous fragments to the atom with max links
1032
  d=1e9
1033
  DO J=1,izm
1034
    Call vecteur(iat,ind_zmat(j,1),x,y,z,vx1,vy1,vz1,norm1)
1035
    if (norm1.le.d) THEN
1036
      Jat=j
1037
      d=norm1
1038
    END IF
1039
  END DO
1040
  n2=ind_zmat(jat,1)
1041
  SELECT CASE(jat)
1042
        CASE (1)
1043
     n3=ind_zmat(2,1)
1044
     n4=ind_zmat(3,1)
1045
    CASE (2)
1046
     n3=ind_zmat(1,1)
1047
     n4=ind_zmat(3,1)
1048
    CASE DEFAULT
1049
     n3=ind_zmat(jAt,2)
1050
     n4=ind_zmat(jat,3)
1051
  END SELECT
1052
  izm=izm+1
1053
  Call add2indzmat(na,izm,iat,n2,n3,n4,ind_zmat,x,y,z)
1054
    idx_zmat(iat)=izm
1055
    DejaFait(iat)=.TRUE.
1056
  IdxCaFaire=2
1057
    CaFaire(1)=iat
1058
              CaFaire(2)=0
1059
              FCaf(Iat)=.TRUE.
1060
              IaFaire=1
1061
              DO WHILE (CaFaire(IaFaire).NE.0)
1062
                 n1=CaFaire(IaFaire)
1063
                 n2=ind_zmat(idx_zmat(N1),2)
1064
                 if (idx_zmat(N1).EQ.2) THEN
1065
                    !     We have a (small) problem: we have to add atoms linked to the 2nd
1066
                    !     atom of the zmat. This is a pb because we do not know
1067
                    !     which atom to use to define the dihedral angle
1068
                    !     we take the third atom of the zmat
1069
                    n3=ind_zmat(3,1)
1070
                 ELSE
1071
                    n3=ind_zmat(idx_zmat(n1),3)
1072
                 END IF
1073
                 DO I3=1,Liaisons(n1,0)
1074
                    IAt=Liaisons(n1,I3)
1075
! PFL 29.Aug.2008 ->
1076
! We dissociate the test on 'DejaFait' that indicates that this atom
1077
! has already been put in the Zmat
1078
! from the test on FCaf that indicates that this atom has been put in the
1079
! 'CAFaire' list that deals with identifying its connectivity.
1080
! Those two test might differ for a frozen atom linked to non frozen atoms.
1081
                    IF (.NOT.DejaFait(Iat)) THEN
1082
                       izm=izm+1
1083
                       Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
1084
                       idx_zmat(iat)=izm
1085
                       n3=iat
1086
                       DejaFait(Iat)=.TRUE.
1087
                    END IF
1088
                    IF (.NOT.FCaf(Iat)) THEN
1089
                       CaFaire(IdxCaFaire)=iat
1090
                       IdxCaFaire=IdxCaFaire+1
1091
                       CaFaire(IdxCaFaire)=0
1092
                       FCaf(Iat)=.TRUE.
1093
                    END IF
1094
! <- PFL 29.Aug.2008
1095
                 END DO
1096
                 IaFaire=IaFaire+1
1097
              END Do              ! DO WHILE CaFaire
1098

    
1099
        if (debug) THEN
1100
           WRITE(*,*) 'ind_zmat 4'
1101
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
1102
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
1103
           WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
1104
           DO Ip=4,izm
1105
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1), ind_zmat(Ip,2), &
1106
                   ind_zmat(Ip,3), ind_zmat(Ip,4)
1107
           END DO
1108
        END IF
1109

    
1110
  END DO                    ! Loop on all fragments
1111

    
1112

    
1113
     ! We have ind_zmat. We calculate val_zmat :-)
1114
     if (debug) WRITE(*,*) "Calculating val_zmat"
1115

    
1116
     val_zmat(1,1)=0.d0
1117
     val_zmat(1,2)=0.d0
1118
     val_zmat(1,3)=0.d0
1119
     val_zmat(2,2)=0.d0
1120
     val_zmat(2,3)=0.d0
1121
     val_zmat(3,3)=0.d0
1122

    
1123
     n1=ind_zmat(2,1)
1124
     n2=ind_zmat(2,2)
1125

    
1126
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
1127

    
1128
     val_zmat(2,1)=norm1
1129

    
1130

    
1131
     n1=ind_zmat(3,1)
1132
     n2=ind_zmat(3,2)
1133
     n3=ind_zmat(3,3)
1134

    
1135
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
1136

    
1137
     val_zmat(3,1)=norm1
1138

    
1139
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
1140
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
1141

    
1142
     val_zmat(3,2)=val
1143

    
1144
     DO i=4,na
1145

    
1146
        n1=ind_zmat(i,1)
1147
        n2=ind_zmat(i,2)
1148
        n3=ind_zmat(i,3)
1149
        n4=ind_zmat(i,4)
1150

    
1151
        if (debug) WRITE(*,*) "Doing i,n1,n2,n3,n4",i,n1,n2,n3,n4
1152
        CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
1153

    
1154
        CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
1155
        val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
1156

    
1157
        CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
1158
        CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2,vx4,vy4,vz4,norm4)
1159
        CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2,vx5,vy5,vz5,norm5)
1160

    
1161
        val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &
1162
             vx2,vy2,vz2,norm2)
1163

    
1164
        !     write(*,11) n1,n2,norm1,n3,val,n4,val_d
1165
        !11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
1166

    
1167
        val_zmat(i,1)=norm1
1168
        val_zmat(i,2)=val
1169
        val_zmat(i,3)=val_d
1170

    
1171
     END DO
1172

    
1173
     if (debug) THEN
1174
        WRITE(*,*) 'DBG Cre_Zmat_Frag: ind_zmat'
1175
        DO I=1,na
1176
           WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
1177
        END DO
1178

    
1179
        WRITE(*,*) 'DBG Cre_Zmat_Frag: Full zmat'
1180
        DO I=1,na
1181
           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)
1182
        END DO
1183

    
1184
     END IF
1185

    
1186
     if (debugGaussian) THEN
1187
        WRITE(*,*) 'DBG Cre_Zmat_Frag: Gaussian Zmat - START'
1188
        Call zmat_g92(na,atome,ind_zmat,val_zmat)
1189
        WRITE(*,*) 'DBG Cre_Zmat_Frag: Gaussian Zmat - END'
1190
     END IF
1191

    
1192

    
1193
     if (debug) WRITE(*,*) "Deallocate (FragDist,Fragment, NbAtFrag,FragAt)"
1194
     DEALLOCATE(FragDist,Fragment, NbAtFrag,FragAt)
1195
     if (debug) WRITE(*,*) "Deallocate (DistFrag,Liaisons)"
1196
     DEALLOCATE(DistFrag,Liaisons)
1197
     if (debug) WRITE(*,*) "Deallocate(CaFaire,DejaFait)"
1198
     DEALLOCATE(CaFaire,DejaFait,FCaf,MaxLFrag)
1199
 
1200

    
1201

    
1202
     if (debug) WRITE(*,*) "=============================== Exiting Calc_zmat_frag ========================"
1203

    
1204
END SUBROUTINE Calc_Zmat_frag
1205

    
1206
SUBROUTINE zmat_g92(na,atome,ind_zmat,val_zmat)
1207

    
1208
! This subroutine comes for Cart. Slightly modified to be f90
1209

    
1210
  Use Path_module, only :  Nom
1211
  Use Io_module
1212

    
1213
  IMPLICIT NONE
1214

    
1215
  integer(KINT), INTENT(IN) ::  na,atome(na)
1216
  INTEGER(KINT), INTENT(IN) :: ind_zmat(Na,5)
1217
  real(KREAL), INTENT(IN)   :: val_zmat(Na,3)
1218

    
1219
  character(6) :: at1, at2, at3, at4, d, a, dh
1220
  character(SCHARS), ALLOCATABLE :: tab(:) ! 3*na
1221
  character(LCHARS) :: ligne
1222
  
1223
  INTEGER(KINT) :: i,n1,n2,n3,n4
1224

    
1225
  ALLOCATE(tab(3*na))
1226

    
1227
  DO i=1,na
1228
     IF (i .GE. 1)  THEN
1229
        n1=ind_zmat(i,1)
1230
        write(at1,11) nom(atome(n1)),n1
1231
11      format(a2,i3)
1232
        Call ecris_sb(at1,at1)
1233
        write(ligne,4) at1
1234
     END IF
1235
     IF (i .GE. 2)  THEN
1236
        n2=ind_zmat(i,2)
1237
        write(at2,11) nom(atome(n2)),n2
1238
        Call ecris_sb(At2,at2)
1239
        write(d,11) 'R',i-1
1240
        Call ecris_sb(D,d)
1241
        write(ligne,4) at1,at2,d
1242
        write(tab(i-1),12) d,val_zmat(i,1)
1243
12      format(a8,f8.4)
1244
     END IF
1245
     IF (i .GE. 3)  THEN
1246
        n3=ind_zmat(i,3)
1247
        write(at3,11) nom(atome(n3)),n3
1248
        Call ecris_sb(At3,at3)
1249
        write(a,11) 'A',na+i-3
1250
        Call ecris_sb(A,A)
1251
        write(ligne,4) at1,at2,d,at3,a
1252
           write(tab(na+i-3),13) a,val_zmat(i,2)
1253
13         format(a8,f8.3)
1254
        END IF
1255
        IF (i .GE. 4)  THEN
1256
           n4=ind_zmat(i,4)
1257
           write(at4,11) nom(atome(n4)),n4
1258
           Call ecris_sb(At4,at4)
1259
           write(dh,11) 'DH',na+na+i-6
1260
           Call ecris_sb(dh,dh)
1261
           write(ligne,4) at1,at2,d,at3,a,at4,dh
1262
4          format(7a6)
1263
           write(tab(na+na+i-6),13) dh,val_zmat(i,3)
1264
        END IF
1265
        write(IOOUT,*) TRIM(ligne)
1266
     END DO
1267
     
1268
     write(IOOUT,*)
1269
     IF (na .EQ. 2) THEN
1270
        write(IOOUT,*) TRIM(tab(1))
1271
     ELSE
1272
        DO i=1,3*na-6
1273
           write(IOOUT,*) TRIM(tab(i))
1274
        END DO
1275
     END IF
1276
     write(IOOUT,*)
1277

    
1278
     DEALLOCATE(Tab)
1279

    
1280
   END SUBROUTINE zmat_g92
1281

    
1282
   SUBROUTINE ecris_sb(inter1,inter)
1283

    
1284
      character(6) :: inter,inter1
1285

    
1286

    
1287
      k=1
1288
      DO j=1,len(inter)
1289
        IF (inter(j:j) .NE. ' ' ) THEN
1290
         inter1(k:k)=inter(j:j)
1291
         k=k+1
1292
        END IF
1293
      END DO
1294

    
1295
       inter1(k:6)=' '
1296

    
1297
    END SUBROUTINE ecris_sb