Statistiques
| Révision :

root / src / Calc_mixed_frag.f90 @ 8

Historique | Voir | Annoter | Télécharger (37,67 ko)

1
SUBROUTINE Calc_mixed_frag(na,atome,ind_zmat,val_zmat,x,y,z, &
2
     r_cov,fact,frozen,cart,ncart)
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
! Parameters of the subroutine
18
! na: number of atoms in the system
19
  integer(KINT) ::  na
20
! atome: Mass number of the atoms of the system
21
  integer(KINT) ::  atome(na)
22
! ind_zmat: for "zmat" atoms contains the indices of reference atoms
23
  integer(KINT) ::  ind_zmat(Na,5)
24
  INTEGER(KINT) ::  idx_zmat(NA)
25
  real(KREAL) ::  x(Na),y(Na),z(Na),fact
26
  real(KREAL) ::  val_zmat(Na,3)
27
  real(KREAL) ::  r_cov(0:Max_Z)
28

    
29
  CHARACTER(5) :: AtName
30

    
31
  !     Frozen contains the indices of frozen atoms
32
  INTEGER(KINT) Frozen(*),Cart(*),NFroz,NCart
33
  LOGICAL, ALLOCATABLE :: FCart(:) ! Na
34
  INTEGER(KINT), ALLOCATABLE :: AtCart(:) !Na
35
  INTEGER(KINT), ALLOCATABLE :: FrozDist(:) ! na
36
  INTEGER(KINT) NbFrag,IdxFrag
37
  INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !na
38
  INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !na,na
39
  INTEGER(KINT), ALLOCATABLE :: FrozFrag(:,:) !na,3
40
  INTEGER(KINT), ALLOCATABLE :: FrozBlock(:,:) !(na,0:na)
41
  REAL(KREAL), ALLOCATABLE :: DistFroz(:) !na
42

    
43
  INTEGER(KINT) :: IdxCaFaire, IAfaire
44
  INTEGER(KINT), ALLOCATABLE ::  LIAISONS(:,:) ! (Na,0:NMaxL)
45
  INTEGER(KINT), ALLOCATABLE ::  LiaisonsBis(:,:) ! (Na,0:NMaxL)
46
  INTEGER(KINT), ALLOCATABLE ::  LiaisonsIni(:,:) ! (Na,0:NMaxL)
47
  INTEGER(KINT), ALLOCATABLE :: CAFaire(:) ! (Na)
48

    
49

    
50
  real(KREAL) ::  vx1,vy1,vz1,norm1
51
  real(KREAL) ::  vx2,vy2,vz2,norm2
52
  real(KREAL) ::  vx3,vy3,vz3,norm3
53
  real(KREAL) ::  vx4,vy4,vz4,norm4
54
  real(KREAL) ::  vx5,vy5,vz5,norm5
55
  real(KREAL) ::  val, val_d, d
56
  Logical Debug,DebugGaussian
57
  LOGICAL, ALLOCATABLE :: DejaFait(:),FCaf(:) !(na)
58
  Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen
59

    
60

    
61

    
62
  INTEGER(KINT) :: I, J, n1, n2, n3, n4, IAt, IL, JL, IFrag, ITmp
63
  INTEGER(KINT) :: I3, I1, Ip, IFragFroz, IBeg
64
  INTEGER(KINT) :: I0, Izm, JAt,Izm0,Idx
65

    
66
  REAL(KREAL) :: Pi, Sang
67

    
68
  INTERFACE
69
     function valid(string) result (isValid)
70
       CHARACTER(*), intent(in) :: string
71
       logical                  :: isValid
72
     END function VALID
73

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

    
81
     FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
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) ::  SinAngle
86
     END FUNCTION SINANGLE
87

    
88

    
89
     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
90
       use Path_module, only :  Pi,KINT, KREAL
91
       real(KREAL) ::  v1x,v1y,v1z,norm1
92
       real(KREAL) ::  v2x,v2y,v2z,norm2
93
       real(KREAL) ::  v3x,v3y,v3z,norm3 
94
       real(KREAL) ::  angle_d,ca,sa
95
     END FUNCTION ANGLE_D
96

    
97
  END INTERFACE
98

    
99
  Pi=dacos(-1.d0)
100
  debug=valid("calc_mixed_frag")
101
  debugGaussian=valid("zmat_gaussian")
102

    
103
  if (debug) Call Header("Entering Calc_mixed_frag")
104
  if (na.le.2) THEN
105
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
106
     RETURN
107
  END IF
108

    
109
  ALLOCATE(FrozDist(Na),Fragment(na), NbAtFrag(na),FragAt(na,na))
110
  ALLOCATE(FCart(Na),AtCart(Na))
111
  ALLOCATE(FrozFrag(na,3), FrozBlock(na,0:na))
112
  ALLOCATE(DistFroz(na),Liaisons(na,0:NMaxL))
113
  ALLOCATE(LiaisonsBis(na,0:NMaxL),LiaisonsIni(na,0:NMaxL))
114
  ALLOCATE(CaFaire(na+1),DejaFait(Na),FCaf(na),FrozAt(na))
115

    
116
  Ind_Zmat=0
117

    
118
  ! Putting cart+frozen atoms into  cart array
119
  FCart=.FALSE.
120
  NCart=0
121
  Idx=1
122
  DO WHILE (Frozen(Idx).NE.0)
123
     If (Frozen(Idx).LT.0) THEN
124
        DO I=Frozen(Idx-1),abs(Frozen(Idx))
125
           IF (.NOT.FCart(I)) THEN
126
              NCart=NCart+1
127
              AtCart(NCart)=I
128
              FCart(I)=.TRUE.
129
           END IF
130
        END DO
131
     ELSEIF (.NOT.FCart(Frozen(Idx))) THEN
132
        FCart(Frozen(Idx))=.TRUE.
133
        NCart=NCart+1
134
        AtCart(NCart)=Frozen(Idx)
135
     END IF
136
  Idx=Idx+1
137
  END DO
138
  NFroz=NCart
139
  Idx=1
140
  DO WHILE (Cart(Idx).NE.0)
141
     If (Cart(Idx).LT.0) THEN
142
        DO I=Cart(Idx-1),abs(Cart(Idx))
143
           IF (.NOT.FCart(I)) THEN
144
              NCart=NCart+1
145
              AtCart(NCart)=I
146
              FCart(I)=.TRUE.
147
           END IF
148
        END DO
149
     ELSEIF (.NOT.FCart(Cart(Idx))) THEN
150
        FCart(Cart(Idx))=.TRUE.
151
        NCart=NCart+1
152
        AtCart(NCart)=Cart(Idx)
153
     END IF
154
  Idx=Idx+1
155
  END DO
156

    
157
  if (debug) THEN
158
     WRITE(*,'(1X,A,I4,A,I4,A)') "Found ",NFroz," frozen atoms, and a total of ",NCart," atoms described in cartesian coordinates"
159
     WRITE(*,*) "AtCart:",AtCart(1:NCart)
160
  END IF
161

    
162
  DejaFait=.FALSE.
163

    
164
  ! We cheat a lot: we will use ind_zmat et val_zmat to store
165
  ! the cartesian coordinates of the NCart atoms :-p
166
  DO I=1,NCart
167
     Iat=AtCart(I)
168
     ind_zmat(I,1)=Iat
169
     ind_zmat(I,2:5)=-1
170
     val_zmat(I,1)=X(Iat)
171
     val_zmat(I,2)=Y(Iat)
172
     val_zmat(I,3)=Z(Iat)
173
     DejaFait(Iat)=.TRUE.
174
     idx_zmat(Iat)=I
175
  END DO
176

    
177

    
178
  izm=NCart
179

    
180
  ! We now calculate the connections
181
  Liaisons=0
182
  LiaisonsBis=0
183

    
184
  if (debug) THEN
185
     WRITE(*,*) "Liaisons initialized"
186
!     WRITE(*,*) 'Covalent radii used'
187
!     DO iat=1,na
188
!        i=atome(iat)
189
!        WRITE(*,*) Nom(I),Iat,r_cov(i),r_cov(i)*fact
190
!     END DO
191
  END IF
192

    
193
1003 FORMAT(1X,I4,' - ',25(I5))
194

    
195
  Call CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)
196

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

    
204
  ! We get rid of bonds between cart atoms
205
  ! the trick is to do this on liaisons while keeping LiaisonsBis ok.
206

    
207
  LiaisonsIni=Liaisons
208
  LiaisonsBis=Liaisons
209

    
210
  DO I=1,NCart
211
     Iat=AtCart(I)
212
     if (debug) WRITE(*,'(20(1X,I3))') I,Iat,  &
213
          (LiaisonsIni(Iat,Itmp),ITmp=0,LiaisonsIni(Iat,0))
214
     DO J=1,LiaisonsIni(Iat,0)
215
        Jat=LiaisonsIni(Iat,J)
216
        IF (FCart(Jat)) THEN
217
           Call Annul(Liaisons,Iat,Jat)
218
           Call Annul(Liaisons,Jat,Iat)
219
           Call Annul(LiaisonsIni,Jat,Iat)
220
           Liaisons(Iat,0)=Liaisons(Iat,0)-1
221
           Liaisons(Jat,0)=Liaisons(Jat,0)-1
222
           LiaisonsIni(Jat,0)=LiaisonsIni(Jat,0)-1
223
        END IF
224
     END DO
225
  END DO
226

    
227
  if (debug) THEN
228
     WRITE(*,*) "Connections omitting bonds between cart atoms"
229
     DO IL=1,na
230
        WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
231
     END DO
232
  END IF
233

    
234
  ! We decompose the system into fragments, but we omit on purpuse
235
  ! fragments consisting only of frozen atoms
236
  ! Now, frozblock will contain the frozen atom indices in a given fragment
237
  Fragment=0
238
  NbAtFrag=0
239
  FrozBlock=0
240

    
241
  IdxFrag=0
242
  NbFrag=0
243

    
244
 Call Decomp_frag(na,liaisons,.NOT.Fcart,nbfrag,Fragment,NbAtFrag,FragAt)
245

    
246
 ! We compute FrozBlock now
247
  DO IFrag=1,NbFrag
248
     DO I=1,NbAtFrag(IFrag)
249
        Iat=FragAt(IFrag,I)
250
        IF (FCart(Iat)) THEN
251
           FrozBlock(IFrag,0)=FrozBlock(IFrag,0)+1
252
           FrozBlock(IFrag,FrozBlock(IFrag,0))=IAt
253
        END IF
254
     END DO
255
  END DO
256

    
257
  if (debug) THEN
258
     WRITE(*,*) 'I found ',NbFrag, 'fragments with respectively '
259
     WRITE(*,*) (NbAtFrag(I),I=1,NbFrag), 'atoms'
260
     WRITE(*,*) "Fragments atoms are now listed as F in the following"
261
     DO I=1,NbFrag
262
        WRITE(*,*) Na
263
        WRITE(*,*) 'Fragment ', i
264
        DO J=1,Na
265
           AtName=Nom(Atome(J))
266
           IF (Fragment(J).EQ.I) AtName='F'
267
           WRITE(*,'(1X,A5,3(1X,F10.6))') AtName,X(J),Y(J),Z(J)
268
        END DO
269
        !     WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
270
     END DO
271

    
272
     DO I=1,NbFrag
273
        WRITE(*,*) 'Fragment ', i
274
        WRITE(*,'(A,20(1X,I3))') "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
275
        WRITE(*,'(A,20(1X,I3))') "FrozBlock:",(FrozBlock(I,j),j=0,FrozBlock(I,0))
276
     END DO
277
  END IF
278

    
279
  NFroz=0
280
  DO IFrag=1,NbFrag
281
     If (FrozBlock(IFrag,0).NE.0) NFroz=NFroz+1
282
  END DO
283

    
284
  IF (DEBUG) WRITE(*,'(1X,A,I4,A,I4,A)') "Found ",NFroz," fragment(s) with cart atoms, out of",NbFrag," fragment(s)"
285

    
286

    
287
  if (debug) THEN
288
     WRITE(*,*) "Connections before adding fragments"
289
     DO IL=1,na
290
        IF (Liaisons(IL,0).NE.0)  WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
291
     END DO
292
  END IF
293

    
294
  ! We will now add the fragments containing non cart atoms.
295
  ! I am not sure that there is a best strategy to add them...
296
  ! so we add them in the order they were created :-(
297
  ! First only block with frozen atoms are added
298
  izm0=-1
299
  IFrag=1
300
  IZm=NCart
301
  FCaf=.FALSE.
302
  DO IFragFroz=1,NFroz
303
     DO WHILE (FrozBlock(IFrag,0).EQ.0)
304
        IFrag=IFrag+1
305
     END DO
306
     ! each atom has at least one frozen atom that will serve as the seed
307
     ! of this fragment.
308
     if (debug) WRITE(*,*) 'Adding fragment :', ifrag,'with',FrozBlock(IFrag,0),' frozen atoms'
309
     IF (FrozBlock(IFrag,0).EQ.1) THEN
310
        ! There is only one frozen atom, easy case ...
311
        Iat=FrozBlock(IFrag,1)
312
        if (debug) WRITE(*,*) 'Only frozen atom of frag',iat
313
        DejaFait(iat)=.TRUE.
314
        IdxCaFaire=2
315
        CaFaire(1)=iat
316
        CaFaire(2)=0
317
        IaFaire=1
318
        FCaf(Iat)=.TRUE.
319
        DO WHILE (CaFaire(IaFaire).NE.0)
320
           n1=CaFaire(IaFaire)
321
           I1=idx_zmat(n1)
322
           n2=ind_zmat(I1,2)
323
           n3=ind_zmat(I1,3)
324

    
325
        if (debug) WRITE(*,1003) n1,(LIAISONS(n1,JL),JL=0,NMaxL)
326
           DO I3=1,Liaisons(n1,0)
327
              IAt=Liaisons(n1,I3)
328
              ! PFL 2007.Apr.16 ->
329
              ! We add ALL atoms linked to n1 to CaFaire
330
              ! But we delete their link to n1
331
!! PFL 2007.Aug.28 ->
332
! re-add the test on FCaf in order not to put the same atom more than once in 
333
! CaFaire array.
334
                 IF (.NOT.FCaf(Iat)) THEN
335
                    CaFaire(IdxCaFaire)=Iat
336
                    IdxCaFaire=IdxCaFaire+1
337
                    CaFaire(IdxCaFaire)=0
338
                    FCaf(Iat)=.TRUE.
339
                 END IF
340
!! <-- PFL 2007.Aug.28
341
              if (debug) WRITE(*,1003) Iat,(LIAISONS(Iat,JL),JL=0,NMaxL)
342
              Call Annul(Liaisons,Iat,n1)
343
              if (debug) WRITE(*,1003) Iat,(LIAISONS(Iat,JL),JL=0,NMaxL)
344
              Liaisons(iat,0)=Liaisons(Iat,0)-1
345
           if (debug) WRITE(*,*) 'Iafaire,n1,DejaFait,idxzmat,n2,n3',IaFaire,n1, &
346
                DejaFait(n1),idx_zmat(n1),n2,3
347
           if (debug) WRITE(*,*) 'Cafaire - 0', &
348
                (CaFaire(J),J=Iafaire,IdxCAfaire)
349

    
350

    
351
              ! <- PFL 2007.Apr.16
352
              IF (.NOT.DejaFait(Iat)) THEN
353
                 ! we add it to the Zmat
354
!                 WRITE(*,*) "coucou"
355
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
356
!                 WRITE(*,*) "coucou 2"
357
                 izm=izm+1
358
                 IF (izm.EQ.2) THEN
359
!                 WRITE(*,*) "coucou 3"
360
                    ind_zmat(izm,1)=IAt
361
                    ind_zmat(izm,2)=n1
362
                    ind_zmat(izm,3:5)=0
363
                    if (n2.EQ.-1) n2=Iat
364
                 END IF
365

    
366
!                 WRITE(*,*) "coucou 4"
367
                 if ((izm.GE.3).AND.(n2.EQ.-1)) THEN
368
!                 WRITE(*,*) "coucou 5"
369
!                 WRITE(*,*) "coucou 3 bis"
370
                       ! This atom is defined in cartesian coordinates
371
                       ! and has not been used to put a non cart atom.
372
                       ! we will find the closest atom linked to this one amongst the atoms
373
                       ! already stored in ind_zmat
374
                       d=1e9
375
                       JAt=-1
376
                       DO I=1,NCart
377
                          If (ind_zmat(I,1).NE.n1) THEN
378
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
379
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
380
                             ! we take only atoms that are not aligned
381
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
382
                                Jat=ind_zmat(I,1)
383
                                d=norm2
384
                             END IF
385
                          END IF
386
                       END DO
387
                       if (JAt.EQ.-1) THEN
388
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
389
                          STOP
390
                       END IF
391
                       n2=JAt
392
                 END IF! izm>=3 and n2.eq.-1
393

    
394
                 IF (izm.eq.3) THEN
395
                    ind_zmat(izm,1)=Iat
396
                    ind_zmat(izm,2)=n1
397
                    ind_zmat(izm,3)=n2
398
                 END IF  ! izm=3
399

    
400
                 IF (izm.GE.4) THEN
401
!                                     WRITE(*,*) "coucou 5 bis"
402
                    IF (n3.EQ.-1) THEN
403
                       ! This atom is defined in cartesian coordinates
404
                       ! and has not been used to put a non cart atom.
405
                       ! we will find the closest atom linked to this one amongst the atoms
406
                       ! already stored in ind_zmat
407
                       call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
408
                       d=1e9
409
                       JAt=-1
410
                       DO I=1,NCart
411
                          If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
412
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
413
                             SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
414
                             if (debug) WRITE(*,*) "1-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
415
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
416
                             ! we take only atoms that are not aligned
417
                             if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
418
                                Jat=ind_zmat(I,1)
419
                                d=norm3
420
                             END IF
421
                          END IF
422
                       END DO
423
                       if (JAt.EQ.-1) THEN
424
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
425
                          STOP
426
                       END IF
427
                       n3=JAt
428
                    END IF
429
                    !                    ind_zmat(I1,3)=n3
430
                    Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
431
                 END IF ! izm>=4
432
!                 WRITE(*,*) "coucou 6"
433
           if (debug) THEN
434
              WRITE(*,*) 'ind_zmat 0',izm
435
              DO Ip=1,NCart
436
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(ip,1:4)
437
                 END DO
438
              SELECT CASE (NCart)
439
                 CASE (1)
440
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
441
                    if (izm.GE.3) &
442
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
443
                    Ibeg=4
444
                 CASE (2)
445
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
446
                    Ibeg=4
447
                 CASE DEFAULT
448
                    IBeg=NCart+1
449
                 END SELECT
450
              DO Ip=IBeg,izm
451
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
452
                      ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
453
              END DO
454
           END IF
455

    
456
                 idx_zmat(iat)=izm
457
!                 Call Annul(Liaisons,iat,n1)
458
!                 Liaisons(iat,0)=Liaisons(Iat,0)-1
459
                 !                  Call Annul(Liaisons,n1,iat)
460
                 n3=iat
461
                 DejaFait(Iat)=.TRUE.
462
              END IF
463
           END DO
464
           IaFaire=IaFaire+1
465

    
466
           if (debug) THEN
467
              WRITE(*,*) 'ind_zmat 5'
468
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
469
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
470
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
471
              DO Ip=4,izm
472
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
473
                      ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
474
              END DO
475
           END IF
476

    
477
        END Do              ! DO WHILE CaFaire
478

    
479

    
480
     ELSE     ! FrozBlock(I,0)==1
481
        if (debug) WRITE(*,*) 'IFrag=',IFrag,'Frozblock(I,0)>1',Frozblock(IFrag,0)
482
        ! We have many frozen atoms for one fragment. We will have to choose
483
        ! the first one, and then to decide when to swich...
484
        Iat=0
485
        I0=-1
486
        DO I=1,FrozBlock(IFrag,0)
487
           JAt=FrozBlock(IFrag,I)
488
           if (debug) WRITE(*,*) Jat, &
489
                (LiaisonsBis(Jat,Itmp),Itmp=0,LiaisonsBis(Jat,0))
490
           IF (LiaisonsBis(Jat,0).GT.I0) THEN
491
              I0=LiaisonsBis(Jat,0)
492
              Iat=Jat
493
           END IF
494
        END DO
495
        ! Iat is the starting frozen atom
496
        IF (debug) WRITE(*,*) 'Choosing atom ',iat,'as a starting atom for frag',ifrag
497
        DejaFait(iat)=.TRUE.
498
        IdxCaFaire=2
499
        CaFaire(1)=iat
500
        CaFaire(2)=0
501
        IaFaire=1
502
        FCaf(Iat)=.TRUE.
503
        DO WHILE (CaFaire(IaFaire).NE.0)
504
           n1=CaFaire(IaFaire)
505
           I1=idx_zmat(n1)
506
           n2=ind_zmat(I1,2)
507
           n3=ind_zmat(I1,3)
508
           if (debug) WRITE(*,*) 'Iafaire,n1,DejaFait,idxzmat,n2,n3',IaFaire,n1, &
509
                DejaFait(n1),idx_zmat(n1),n2,3
510
           if (debug) WRITE(*,*) 'Cafaire', &
511
                (CaFaire(J),J=Iafaire,IdxCAfaire)
512

    
513

    
514
           DO I3=1,Liaisons(n1,0)
515
              if (debug) WRITE(*,*) "Here, I3=",I3
516
              IAt=Liaisons(n1,I3)
517
              ! PFL 2007.Apr.16 ->
518
              ! We add ALL atoms linked to n1 to CaFaire
519
              ! But we delete their link to n1
520
!! PFL 2007.Aug.28 ->
521
! re-add the test on FCaf in order not to put the same atom more than once in 
522
! CaFaire array.
523
                 IF (.NOT.FCaf(Iat)) THEN
524
                    CaFaire(IdxCaFaire)=Iat
525
                    IdxCaFaire=IdxCaFaire+1
526
                    CaFaire(IdxCaFaire)=0
527
                    FCaf(Iat)=.TRUE.
528
                 END IF
529
!! <-- PFL 2007.Aug.28
530
              Call Annul(Liaisons,Iat,n1)
531
              Liaisons(iat,0)=Liaisons(Iat,0)-1
532
              ! <- PFL 2007.Apr.16
533
              IF (.NOT.DejaFait(Iat)) THEN
534
                 ! we add it to the Zmat
535
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
536
                 izm=izm+1
537
                 IF (izm.EQ.2) THEN
538
                    ind_zmat(izm,1)=IAt
539
                    ind_zmat(izm,2)=n1
540
                    ind_zmat(izm,3:5)=0
541
                 ELSE
542
                    IF (n2.EQ.-1) THEN
543
                       ! This atom is defined in cartesian coordinates
544
                       ! and has not been used to put a non cart atom.
545
                       ! we will find the closest atom linked to this one amongst the atoms
546
                       ! already stored in ind_zmat
547
                       d=1e9
548
                       JAt=-1
549
                       DO I=1,NCart
550
                          If (ind_zmat(I,1).NE.n1) THEN
551
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
552
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
553
                             ! we take only atoms that are not aligned
554
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
555
                                Jat=ind_zmat(I,1)
556
                                d=norm2
557
                             END IF
558
                          END IF
559
                       END DO
560
                       if (JAt.EQ.-1) THEN
561
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
562
                          STOP
563
                       END IF
564
                       n2=JAt
565
                    END IF
566
                    !                    ind_zmat(I1,2)=n2
567
                    if (izm.EQ.3) THEN
568
                       ind_zmat(izm,1)=Iat
569
                       ind_zmat(izm,2)=n1
570
                       ind_zmat(izm,3)=n2
571
                    ELSE
572
                       IF (n3.EQ.-1) THEN
573
                          ! This atom is defined in cartesian coordinates
574
                          ! and has not been used to put a non cart atom.
575
                          ! we will find the closest atom linked to this one amongst the atoms
576
                          ! already stored in ind_zmat
577
                          call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
578
                          d=1e9
579
                          JAt=-1
580
                          DO I=1,NCart
581
                             If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
582
                                call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
583
                                SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
584
                             if (debug) WRITE(*,*) "2-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
585
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
586

    
587
                                ! we take only atoms that are not aligned
588
                                if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
589
                                   Jat=ind_zmat(I,1)
590
                                   d=norm3
591
                                END IF
592
                             END IF
593
                          END DO
594
                          if (JAt.EQ.-1) THEN
595
                             WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
596
                             STOP
597
                          END IF
598
                          n3=JAt
599
                       END IF
600
                       !                    ind_zmat(I1,3)=n3
601
                       Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
602
                    END IF
603
                 END IF
604
                 idx_zmat(iat)=izm
605
              !                  Call Annul(Liaisons,n1,iat)
606
                 n3=iat
607
                 DejaFait(Iat)=.TRUE.
608
                 if (debug) WRITE(*,*) "For no reason, I3=",I3
609
              END IF
610
           END DO
611
           IaFaire=IaFaire+1
612

    
613
           if (debug) THEN
614
              WRITE(*,*) 'ind_zmat 6',izm
615
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
616
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
617
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
618
              DO Ip=4,izm
619
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
620
                      ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
621
              END DO
622
           END IF
623

    
624

    
625
        END Do              ! DO WHILE CaFaire
626

    
627

    
628
     END IF  ! FrozBlock(I,0)==1
629

    
630
     IFrag=IFrag+1
631

    
632
  END DO                    ! Loop on all fragments
633

    
634

    
635
  DO IFrag=1,NbFrag
636
     IF (FrozBlock(IFrag,0).EQ.0) THEN
637
        if (debug) WRITE(*,*) 'Adding fragment :', ifrag,'with no frozen atoms'
638
        ! We have no frozen atoms for this fragment. We will have to choose
639
        ! the first atom to put.
640
        ! For now, we do not care of the 'central' atom (ie the one with
641
        ! no CP atoms...)
642
        ! We just take the atom that is the closest to the already placed
643
        ! atoms !
644

    
645
        I0=0
646
        I1=0
647
        d=1e9
648
        DO J=1,izm
649
           Jat=ind_zmat(J,1)
650
           DO I=1,NbAtfrag(IFrag)
651
              IAt=FragAt(IFrag,I)
652
              IF (.NOT.DejaFait(Iat)) THEN
653
                 Call vecteur(Iat,Jat,x,y,z,vx1,vy1,vz1,norm1)
654
                 IF (norm1.LT.d) THEN
655
                    I1=Jat
656
                    I0=Iat
657
                    d=Norm1
658
                 END IF
659
              END IF
660
           END DO
661
        END DO
662

    
663
        Iat=I0
664
        n1=I1
665
        I1=idx_zmat(n1)
666
        n2=ind_zmat(I1,2)
667
        n3=ind_zmat(I1,3)
668

    
669
        IF (debug) WRITE(*,*) Iat,' starting atom for frag',ifrag
670

    
671
        
672
        call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
673
        izm=izm+1
674
        IF (izm.EQ.2) THEN
675
           ind_zmat(izm,1)=IAt
676
           ind_zmat(izm,2)=n1
677
           ind_zmat(izm,3:5)=0
678
        ELSE
679
           IF (n2.EQ.-1) THEN
680
! This atom is defined in cartesian coordinates
681
! and has not been used to put a non cart atom.
682
! we will find the closest atom linked to this one amongst the atoms
683
! already stored in ind_zmat
684
              d=1e9
685
              JAt=-1
686
              DO I=1,NCart
687
                 If (ind_zmat(I,1).NE.n1) THEN
688
                    call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
689
                    SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
690
! we take only atoms that are not aligned
691
                    if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
692
                       Jat=ind_zmat(I,1)
693
                       d=norm2
694
                    END IF
695
                 END IF
696
              END DO
697
              if (JAt.EQ.-1) THEN
698
                 WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
699
                 STOP
700
              END IF
701
              n2=JAt
702
           END IF
703
!                    ind_zmat(I1,2)=n2
704
        END IF
705
        if (izm.EQ.3) THEN
706
           ind_zmat(izm,1)=Iat
707
           ind_zmat(izm,2)=n1
708
           ind_zmat(izm,3)=n2
709
        ELSE
710
           IF (n3.EQ.-1) THEN
711
! This atom is defined in cartesian coordinates
712
! and has not been used to put a non cart atom.
713
! we will find the closest atom linked to this one amongst the atoms
714
! already stored in ind_zmat
715
              call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
716
              d=1e9
717
              JAt=-1
718
              DO I=1,NCart
719
                 If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
720
                    call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
721
                    SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
722
                             if (debug) WRITE(*,*) "3-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
723
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
724

    
725
! we take only atoms that are not aligned
726
                    if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
727
                       Jat=ind_zmat(I,1)
728
                       d=norm3
729
                    END IF
730
                 END IF
731
              END DO
732
              if (JAt.EQ.-1) THEN
733
                 WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
734
                 STOP
735
              END IF
736
              n3=JAt
737
           END IF
738
           !                    ind_zmat(I1,3)=n3
739
           Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
740
        END IF
741
        idx_zmat(iat)=izm
742

    
743

    
744
        DejaFait(iat)=.TRUE.
745
        IdxCaFaire=2
746
        CaFaire(1)=iat
747
        CaFaire(2)=0
748
        IaFaire=1
749
        FCaf(Iat)=.TRUE.
750
        DO WHILE (CaFaire(IaFaire).NE.0)
751
           n1=CaFaire(IaFaire)
752
           if (debug) WRITE(*,*) 'Iafaire,n1',IaFaire,n1, &
753
                DejaFait(n1),idx_zmat(n1)
754
           if (debug) WRITE(*,*) 'Cafaire', &
755
                (CaFaire(J),J=Iafaire,IdxCAfaire)
756
           I1=idx_zmat(n1)
757
           n2=ind_zmat(I1,2)
758
           n3=ind_zmat(I1,3)
759

    
760
           DO I3=1,Liaisons(n1,0)
761
              IAt=Liaisons(n1,I3)
762
              ! PFL 2007.Apr.16 ->
763
              ! We add ALL atoms linked to n1 to CaFaire
764
              ! But we delete their link to n1
765
!! PFL 2007.Aug.28 ->
766
! re-add the test on FCaf in order not to put the same atom more than once in 
767
! CaFaire array.
768
                 IF (.NOT.FCaf(Iat)) THEN
769
                    CaFaire(IdxCaFaire)=Iat
770
                    IdxCaFaire=IdxCaFaire+1
771
                    CaFaire(IdxCaFaire)=0
772
                    FCaf(Iat)=.TRUE.
773
                 END IF
774
!! <-- PFL 2007.Aug.28
775
              Call Annul(Liaisons,Iat,n1)
776
              Liaisons(iat,0)=Liaisons(Iat,0)-1
777
              ! <- PFL 2007.Apr.16
778
              IF (.NOT.DejaFait(Iat)) THEN
779
                 ! we add it to the Zmat
780
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
781
                 izm=izm+1
782
                 IF (izm.EQ.2) THEN
783
                    ind_zmat(izm,1)=IAt
784
                    ind_zmat(izm,2)=n1
785
                    ind_zmat(izm,3:5)=0
786
                 ELSE
787
                    IF (n2.EQ.-1) THEN
788
! This atom is defined in cartesian coordinates
789
! and has not been used to put a non cart atom.
790
! we will find the closest atom linked to this one amongst the atoms
791
! already stored in ind_zmat
792
                       d=1e9
793
                       JAt=-1
794
                       DO I=1,NCart
795
                          If (ind_zmat(I,1).NE.n1) THEN
796
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
797
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
798
! we take only atoms that are not aligned
799
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
800
                                Jat=ind_zmat(I,1)
801
                                d=norm2
802
                             END IF
803
                          END IF
804
                       END DO
805
                       if (JAt.EQ.-1) THEN
806
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
807
                          STOP
808
                       END IF
809
                       n2=JAt
810
                    END IF
811
!                    ind_zmat(I1,2)=n2
812
                 END IF
813
                 if (izm.EQ.3) THEN
814
                    ind_zmat(izm,1)=Iat
815
                    ind_zmat(izm,2)=n1
816
                    ind_zmat(izm,3)=n2
817
                 ELSE
818
                    IF (n3.EQ.-1) THEN
819
! This atom is defined in cartesian coordinates
820
! and has not been used to put a non cart atom.
821
! we will find the closest atom linked to this one amongst the atoms
822
! already stored in ind_zmat
823
                       call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
824
                       d=1e9
825
                       JAt=-1
826
                       DO I=1,NCart
827
                          If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
828
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
829
                             SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
830
                             if (debug) WRITE(*,*) "4-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
831
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
832

    
833
! we take only atoms that are not aligned
834
                             if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
835
                                Jat=ind_zmat(I,1)
836
                                d=norm3
837
                             END IF
838
                          END IF
839
                       END DO
840
                       if (JAt.EQ.-1) THEN
841
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
842
                          STOP
843
                       END IF
844
                       n3=JAt
845
                    END IF
846
!                    ind_zmat(I1,3)=n3
847
                    Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
848
                 END IF
849
                 idx_zmat(iat)=izm
850
                 n3=iat
851
                 DejaFait(Iat)=.TRUE.
852
              END IF
853
           END DO
854
           IaFaire=IaFaire+1
855

    
856
           if (debug) THEN
857
              WRITE(*,*) 'ind_zmat 7',izm
858
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
859
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
860
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
861
              DO Ip=4,izm
862
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
863
                      ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
864
              END DO
865
           END IF
866

    
867

    
868
        END Do              ! DO WHILE CaFaire
869
     END IF                 ! FrozBlock(IFrag,0).EQ.0
870

    
871
  END DO                    ! Loop on all fragments without frozen atoms
872

    
873
  if (debug) THEN
874
     WRITE(*,*) Na+Izm
875
     WRITE(*,*) 'Done... :-(', izm
876
     DO J=1,Na
877
        WRITE(*,'(1X,A5,3(1X,F10.6))') Nom(Atome(J)),X(J),Y(J),Z(J)
878
     END DO
879
     DO I=1,Izm
880
        J=ind_zmat(I,1)
881
        WRITE(*,'(1X,A5,3(1X,F10.6))') 'X ',X(J),Y(J),Z(J)
882
     END DO
883
        IF (izm.NE.Na) THEN
884
           WRITE(*,*) "Atoms not done:"
885
           DO I=1,Na
886
              IF (.NOT.DejaFait(I)) WRITE(*,'(I5)',ADVANCE='NO') I
887
           END DO
888
        END IF
889

    
890
  END IF
891

    
892

    
893
  ! We have ind_zmat. We calculate val_zmat :-)
894
  if (debug) WRITE(*,*) "Calculating val_zmat"
895

    
896
  SELECT CASE (NCart)
897
     CASE (1)
898
        val_zmat(2,2)=0.d0
899
     val_zmat(2,3)=0.d0
900
     val_zmat(3,3)=0.d0
901
     n1=ind_zmat(2,1)
902
     n2=ind_zmat(2,2)
903
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
904
     val_zmat(2,1)=norm1
905

    
906
     n1=ind_zmat(3,1)
907
     n2=ind_zmat(3,2)
908
     n3=ind_zmat(3,3)
909
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
910
     val_zmat(3,1)=norm1
911

    
912
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
913
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
914
     val_zmat(3,2)=val
915
     IBeg=4
916

    
917
     CASE (2)
918
     val_zmat(3,3)=0.d0
919

    
920
     n1=ind_zmat(3,1)
921
     n2=ind_zmat(3,2)
922
     n3=ind_zmat(3,3)
923
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
924
     val_zmat(3,1)=norm1
925

    
926
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
927
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
928
     val_zmat(3,2)=val
929
     IBeg=4
930
  CASE DEFAULT
931
     Ibeg=NCart+1
932
  END SELECT
933

    
934

    
935
  DO i=IBeg,na
936

    
937
     n1=ind_zmat(i,1)
938
     n2=ind_zmat(i,2)
939
     n3=ind_zmat(i,3)
940
     n4=ind_zmat(i,4)
941

    
942
     if (debug) WRITE(*,*) "Doing i,n1,n2,n3,n4",i,n1,n2,n3,n4
943
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
944

    
945
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
946
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
947

    
948
     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
949
     CALL produit_vect(vx1,vy1,vz1,vx2,vy2,vz2, &
950
          vx4,vy4,vz4,norm4)
951
     CALL produit_vect(vx3,vy3,vz3,vx2,vy2,vz2, &
952
          vx5,vy5,vz5,norm5)
953

    
954
     val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &
955
          vx2,vy2,vz2,norm2)
956

    
957
     !     write(*,11) n1,n2,norm1,n3,val,n4,val_d
958
!11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
959

    
960
     val_zmat(i,1)=norm1
961
     val_zmat(i,2)=val
962
     val_zmat(i,3)=val_d
963

    
964
  END DO
965

    
966
  if (debug) THEN
967
     WRITE(*,*) 'DBG Cre_Zmat_Constr: ind_zmat'
968
     DO I=1,na
969
        WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
970
     END DO
971

    
972
     WRITE(*,*) 'DBG Cre_Zmat_Constr: Full zmat'
973
     DO I=1,NCart
974
        WRITE(*,'(1X,I5,3(1X,F11.6))') ind_zmat(i,1),(val_zmat(i,j),j=1,3)
975
     END DO
976
     DO I=NCart+1,Na
977
        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)
978
     END DO
979

    
980
  END IF
981

    
982
  if (debug) WRITE(*,*) "Deallocate (FrozDist,Fragment, NbAtFrag,FragAt)"
983
  DEALLOCATE(FrozDist,Fragment, NbAtFrag,FragAt)
984
  if (debug) WRITE(*,*) "Deallocate FrozFrag, IdxFragAt, FrozBlock"
985
!  DEALLOCATE(FrozFrag, IdxFragAt, FrozBlock)
986
  DEALLOCATE(FrozFrag,FrozBlock)
987
  if (debug) WRITE(*,*) "Deallocate (DistFroz,Liaisons)"
988
  DEALLOCATE(DistFroz,Liaisons)
989
  if (debug) WRITE(*,*) "Deallocate (LiaisonsIni)"
990
  DEALLOCATE(LiaisonsIni)
991
  if (debug) WRITE(*,*) "Deallocate(CaFaire,DejaFait,FrozAt)"
992
  DEALLOCATE(CaFaire,DejaFait,FrozAt)
993
  if (debug) WRITE(*,*) "Deallocate LiaisonsBis"
994
  DEALLOCATE(LiaisonsBis)
995
  if (debug) WRITE(*,*) "Deallocate FCart,AtCart,FCaf"
996
  DEALLOCATE(FCart,AtCart,FCaf)
997

    
998
  if (debugGaussian) THEN
999
     WRITE(*,*) 'DBG Cre_Zmat_Frag: Gaussian Zmat - START'
1000
     Call WriteMixed_Gaussian(na,atome,NCart,ind_zmat,val_zmat)
1001
     WRITE(*,*) 'DBG Cre_Zmat_Frag: Gaussian Zmat - END'
1002
  END IF
1003

    
1004

    
1005
  if (debug) Call Header("Exiting Calc_mixed_frag")
1006

    
1007
END SUBROUTINE Calc_mixed_frag
1008

    
1009

    
1010
SUBROUTINE WriteMixed_Gaussian(na,atome,NCart,ind_zmat,val_zmat)
1011

    
1012
! This subroutine comes for zmat_g92
1013

    
1014
  Use Path_module, only :   Nom
1015
  Use Io_module
1016

    
1017
  IMPLICIT NONE
1018
! Parameters of the subroutine
1019
! na: number of atoms in the system
1020
  integer(KINT), INTENT(IN) ::  na
1021
! atome: Mass number of the atoms of the system
1022
  integer(KINT), INTENT(IN) ::  atome(na)
1023
! nCart: number of atoms described in cartesian
1024
  integer(KINT), INTENT(IN) ::  NCart
1025
! ind_zmat: for "zmat" atoms contains the indices of reference atoms
1026
  integer(KINT), INTENT(IN) ::  ind_zmat(Na,5)
1027
! val_zmat: for "zmat" atoms contains the values from reference atoms
1028
  real(KREAL), INTENT(IN)   :: val_zmat(Na,3)
1029

    
1030
  character(6) :: at1, at2, at3, at4, d, a, dh
1031
  Character(6) :: x,y,z
1032
  character(SCHARS), ALLOCATABLE :: tab(:) ! 3*na
1033
  character(LCHARS) :: ligne
1034
  
1035
  INTEGER(KINT) :: i,n1,n2,n3,n4, NZmat,ITab
1036

    
1037
  ALLOCATE(tab(3*na))
1038

    
1039
! We first write the Cartesian atoms
1040

    
1041
  DO I=1,NCart
1042
! Name of the atom
1043
        n1=ind_zmat(i,1)
1044
        write(at1,11) nom(atome(n1)),n1
1045
        Call ecris_sb(at1,at1)
1046
! x
1047
        write(x,11) 'x',i
1048
        Call ecris_sb(x,x)
1049
        write(tab(3*i-2),'(A,1X,F12.8)') trim(x),val_zmat(i,1)
1050
! y
1051
        write(y,11) 'y',i
1052
        Call ecris_sb(y,y)
1053
        write(tab(3*i-1),'(A,1X,F12.8)') trim(y),val_zmat(i,2)
1054
! z
1055
        write(z,11) 'z',i
1056
        Call ecris_sb(z,z)
1057
        write(tab(3*i),'(A,1X,F12.8)') trim(z),val_zmat(i,3)
1058
! write the line
1059
           write(ligne,'(1X,A,1X,A,1X,A,1X,A)') Trim(at1),x,y,z
1060
        END DO
1061

    
1062
  NZmat=na-NCart
1063
  ITab=3*NCart+1
1064

    
1065
  DO i=NCart+1,NZmat
1066
     IF (i .GE. 1)  THEN
1067
        n1=ind_zmat(i,1)
1068
        write(at1,11) nom(atome(n1)),n1
1069
11      format(a2,i3)
1070
        Call ecris_sb(at1,at1)
1071
        write(ligne,4) at1
1072
     END IF
1073
     IF (i .GE. 2)  THEN
1074
        n2=ind_zmat(i,2)
1075
        write(at2,11) nom(atome(n2)),n2
1076
        Call ecris_sb(At2,at2)
1077
        write(d,11) 'R',i-1
1078
        Call ecris_sb(D,d)
1079
        write(ligne,4) at1,at2,d
1080
        write(tab(ITab),12) d,val_zmat(i,1)
1081
        ITab=ITab+1
1082
12      format(a8,f8.4)
1083
     END IF
1084
     IF (i .GE. 3)  THEN
1085
        n3=ind_zmat(i,3)
1086
        write(at3,11) nom(atome(n3)),n3
1087
        Call ecris_sb(At3,at3)
1088
        write(a,11) 'A',na+i-3
1089
        Call ecris_sb(A,A)
1090
        write(ligne,4) at1,at2,d,at3,a
1091
        write(tab(ITab),13) a,val_zmat(i,2)
1092
        ITab=ITab+1
1093
13      format(a8,f8.3)
1094
     END IF
1095
     IF (i .GE. 4)  THEN
1096
        n4=ind_zmat(i,4)
1097
        write(at4,11) nom(atome(n4)),n4
1098
        Call ecris_sb(At4,at4)
1099
        write(dh,11) 'DH',na+na+i-6
1100
        Call ecris_sb(dh,dh)
1101
        write(ligne,4) at1,at2,d,at3,a,at4,dh
1102
4       format(7a6)
1103
        write(tab(ITab),13) dh,val_zmat(i,3)
1104
        ITab=ITab+1
1105
     END IF
1106
     write(IOOUT,*) TRIM(ligne)
1107
     END DO
1108
     
1109
     write(IOOUT,*)
1110
     DO i=1,ITab
1111
        write(IOOUT,*) TRIM(tab(i))
1112
     END DO
1113
     write(IOOUT,*)
1114

    
1115
     DEALLOCATE(Tab)
1116

    
1117
   END SUBROUTINE WriteMixed_Gaussian