Statistiques
| Révision :

root / src / Calc_mixed_frag.f90 @ 1

Historique | Voir | Annoter | Télécharger (34,46 ko)

1 1 pfleura2
SUBROUTINE Calc_mixed_frag(na,atome,ind_zmat,val_zmat,x,y,z, &
2 1 pfleura2
     r_cov,fact,frozen,cart,ncart)
3 1 pfleura2
4 1 pfleura2
  !
5 1 pfleura2
  !     This second version enables the user to freeze some atoms
6 1 pfleura2
  !     the frozen atoms indices are listed in the frozen array.
7 1 pfleura2
  !
8 1 pfleura2
  !     The idea is surely too stupid to work all the time... but here it is
9 1 pfleura2
  !     we first construct the zmat using only the frozen atoms, and then
10 1 pfleura2
  !     we add the other atoms on top of the first ones...
11 1 pfleura2
  !
12 1 pfleura2
  Use Path_module, only : max_Z, NMaxL, Nom,MaxFroz
13 1 pfleura2
  Use Io_module
14 1 pfleura2
15 1 pfleura2
  IMPLICIT NONE
16 1 pfleura2
17 1 pfleura2
  CHARACTER(5) :: AtName
18 1 pfleura2
  integer(KINT) ::  na,atome(na),ind_zmat(Na,5)
19 1 pfleura2
  INTEGER(KINT) ::  idx_zmat(NA)
20 1 pfleura2
  real(KREAL) ::  x(Na),y(Na),z(Na),fact
21 1 pfleura2
  real(KREAL) ::  val_zmat(Na,3)
22 1 pfleura2
  real(KREAL) ::  r_cov(0:Max_Z)
23 1 pfleura2
24 1 pfleura2
  !     Frozen contains the indices of frozen atoms
25 1 pfleura2
  INTEGER(KINT) Frozen(*),Cart(*),NFroz,NCart
26 1 pfleura2
  LOGICAL, ALLOCATABLE :: FCart(:) ! Na
27 1 pfleura2
  INTEGER(KINT), ALLOCATABLE :: AtCart(:) !Na
28 1 pfleura2
  INTEGER(KINT), ALLOCATABLE :: FrozDist(:) ! na
29 1 pfleura2
  INTEGER(KINT) NbFrag,IdxFrag
30 1 pfleura2
  INTEGER(KINT), ALLOCATABLE :: Fragment(:),NbAtFrag(:) !na
31 1 pfleura2
  INTEGER(KINT), ALLOCATABLE :: FragAt(:,:) !na,na
32 1 pfleura2
  INTEGER(KINT), ALLOCATABLE :: FrozFrag(:,:) !na,3
33 1 pfleura2
  !  INTEGER(KINT), ALLOCATABLE :: IdxFragAt(:) !na
34 1 pfleura2
  INTEGER(KINT), ALLOCATABLE :: FrozBlock(:,:) !(na,0:na)
35 1 pfleura2
  REAL(KREAL), ALLOCATABLE :: DistFroz(:) !na
36 1 pfleura2
37 1 pfleura2
  INTEGER(KINT) :: IdxCaFaire, IAfaire
38 1 pfleura2
  INTEGER(KINT), ALLOCATABLE ::  LIAISONS(:,:) ! (Na,0:NMaxL)
39 1 pfleura2
  INTEGER(KINT), ALLOCATABLE ::  LiaisonsBis(:,:) ! (Na,0:NMaxL)
40 1 pfleura2
  INTEGER(KINT), ALLOCATABLE ::  LiaisonsIni(:,:) ! (Na,0:NMaxL)
41 1 pfleura2
  INTEGER(KINT), ALLOCATABLE :: CAFaire(:) ! (Na)
42 1 pfleura2
43 1 pfleura2
44 1 pfleura2
  real(KREAL) ::  vx1,vy1,vz1,norm1
45 1 pfleura2
  real(KREAL) ::  vx2,vy2,vz2,norm2
46 1 pfleura2
  real(KREAL) ::  vx3,vy3,vz3,norm3
47 1 pfleura2
  real(KREAL) ::  vx4,vy4,vz4,norm4
48 1 pfleura2
  real(KREAL) ::  vx5,vy5,vz5,norm5
49 1 pfleura2
  real(KREAL) ::  val,val_d, d12, d13,d23,d
50 1 pfleura2
  Logical PasFini,Debug
51 1 pfleura2
  LOGICAL, ALLOCATABLE :: DejaFait(:),FCaf(:) !(na)
52 1 pfleura2
  Logical, ALLOCATABLE ::  FrozAt(:)    !T if this atom is frozen
53 1 pfleura2
  LOGICAL F1213, F1223,F1323
54 1 pfleura2
55 1 pfleura2
56 1 pfleura2
  INTEGER(KINT) :: I,J,n0,n1,n2,n3,n4,IAt,IL,JL,IFrag,ITmp, K, KMax
57 1 pfleura2
  INTEGER(KINT) :: IMax, I3,I1,Ip, IFragFroz,IBeg
58 1 pfleura2
  INTEGER(KINT) :: I0, Izm, JAt,Izm0,Idx
59 1 pfleura2
60 1 pfleura2
  REAL(KREAL) :: sDihe, Pi, Sang
61 1 pfleura2
62 1 pfleura2
  INTERFACE
63 1 pfleura2
     function valid(string) result (isValid)
64 1 pfleura2
       CHARACTER(*), intent(in) :: string
65 1 pfleura2
       logical                  :: isValid
66 1 pfleura2
     END function VALID
67 1 pfleura2
68 1 pfleura2
     FUNCTION angle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
69 1 pfleura2
       use Path_module, only :  Pi,KINT, KREAL
70 1 pfleura2
       real(KREAL) ::  v1x,v1y,v1z,norm1
71 1 pfleura2
       real(KREAL) ::  v2x,v2y,v2z,norm2
72 1 pfleura2
       real(KREAL) ::  angle
73 1 pfleura2
     END FUNCTION ANGLE
74 1 pfleura2
75 1 pfleura2
     FUNCTION SinAngle(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2)
76 1 pfleura2
       use Path_module, only :  Pi,KINT, KREAL
77 1 pfleura2
       real(KREAL) ::  v1x,v1y,v1z,norm1
78 1 pfleura2
       real(KREAL) ::  v2x,v2y,v2z,norm2
79 1 pfleura2
       real(KREAL) ::  SinAngle
80 1 pfleura2
     END FUNCTION SINANGLE
81 1 pfleura2
82 1 pfleura2
83 1 pfleura2
     FUNCTION angle_d(v1x,v1y,v1z,norm1,v2x,v2y,v2z,norm2,v3x,v3y,v3z,norm3)
84 1 pfleura2
       use Path_module, only :  Pi,KINT, KREAL
85 1 pfleura2
       real(KREAL) ::  v1x,v1y,v1z,norm1
86 1 pfleura2
       real(KREAL) ::  v2x,v2y,v2z,norm2
87 1 pfleura2
       real(KREAL) ::  v3x,v3y,v3z,norm3
88 1 pfleura2
       real(KREAL) ::  angle_d,ca,sa
89 1 pfleura2
     END FUNCTION ANGLE_D
90 1 pfleura2
91 1 pfleura2
  END INTERFACE
92 1 pfleura2
93 1 pfleura2
  Pi=dacos(-1.d0)
94 1 pfleura2
  debug=valid("calc_mixed_frag")
95 1 pfleura2
  if (debug) Call Header("Entering Calc_mixed_frag")
96 1 pfleura2
  if (na.le.2) THEN
97 1 pfleura2
     WRITE(*,*) "I do not work for less than 2 atoms :-p"
98 1 pfleura2
     RETURN
99 1 pfleura2
  END IF
100 1 pfleura2
101 1 pfleura2
  ALLOCATE(FrozDist(Na),Fragment(na), NbAtFrag(na),FragAt(na,na))
102 1 pfleura2
  ALLOCATE(FCart(Na),AtCart(Na))
103 1 pfleura2
  ALLOCATE(FrozFrag(na,3), FrozBlock(na,0:na))
104 1 pfleura2
  ALLOCATE(DistFroz(na),Liaisons(na,0:NMaxL))
105 1 pfleura2
  ALLOCATE(LiaisonsBis(na,0:NMaxL),LiaisonsIni(na,0:NMaxL))
106 1 pfleura2
  ALLOCATE(CaFaire(na+1),DejaFait(Na),FCaf(na),FrozAt(na))
107 1 pfleura2
108 1 pfleura2
  Ind_Zmat=0
109 1 pfleura2
110 1 pfleura2
  ! Putting cart+frozen atoms into  cart array
111 1 pfleura2
  FCart=.FALSE.
112 1 pfleura2
  NCart=0
113 1 pfleura2
  Idx=1
114 1 pfleura2
  DO WHILE (Frozen(Idx).NE.0)
115 1 pfleura2
     If (Frozen(Idx).LT.0) THEN
116 1 pfleura2
        DO I=Frozen(Idx-1),abs(Frozen(Idx))
117 1 pfleura2
           IF (.NOT.FCart(I)) THEN
118 1 pfleura2
              NCart=NCart+1
119 1 pfleura2
              AtCart(NCart)=I
120 1 pfleura2
              FCart(I)=.TRUE.
121 1 pfleura2
           END IF
122 1 pfleura2
        END DO
123 1 pfleura2
     ELSEIF (.NOT.FCart(Frozen(Idx))) THEN
124 1 pfleura2
        FCart(Frozen(Idx))=.TRUE.
125 1 pfleura2
        NCart=NCart+1
126 1 pfleura2
        AtCart(NCart)=Frozen(Idx)
127 1 pfleura2
     END IF
128 1 pfleura2
  Idx=Idx+1
129 1 pfleura2
  END DO
130 1 pfleura2
  NFroz=NCart
131 1 pfleura2
  Idx=1
132 1 pfleura2
  DO WHILE (Cart(Idx).NE.0)
133 1 pfleura2
     If (Cart(Idx).LT.0) THEN
134 1 pfleura2
        DO I=Cart(Idx-1),abs(Cart(Idx))
135 1 pfleura2
           IF (.NOT.FCart(I)) THEN
136 1 pfleura2
              NCart=NCart+1
137 1 pfleura2
              AtCart(NCart)=I
138 1 pfleura2
              FCart(I)=.TRUE.
139 1 pfleura2
           END IF
140 1 pfleura2
        END DO
141 1 pfleura2
     ELSEIF (.NOT.FCart(Cart(Idx))) THEN
142 1 pfleura2
        FCart(Cart(Idx))=.TRUE.
143 1 pfleura2
        NCart=NCart+1
144 1 pfleura2
        AtCart(NCart)=Cart(Idx)
145 1 pfleura2
     END IF
146 1 pfleura2
  Idx=Idx+1
147 1 pfleura2
  END DO
148 1 pfleura2
149 1 pfleura2
  if (debug) THEN
150 1 pfleura2
     WRITE(*,'(1X,A,I4,A,I4,A)') "Found ",NFroz," frozen atoms, and a total of ",NCart," atoms described in cartesian coordinates"
151 1 pfleura2
     WRITE(*,*) "AtCart:",AtCart(1:NCart)
152 1 pfleura2
  END IF
153 1 pfleura2
154 1 pfleura2
  DejaFait=.FALSE.
155 1 pfleura2
156 1 pfleura2
  ! We cheat a lot: we will use ind_zmat et val_zmat to store
157 1 pfleura2
  ! the cartesian coordinates of the NCart atoms :-p
158 1 pfleura2
  DO I=1,NCart
159 1 pfleura2
     Iat=AtCart(I)
160 1 pfleura2
     ind_zmat(I,1)=Iat
161 1 pfleura2
     ind_zmat(I,2:5)=-1
162 1 pfleura2
     val_zmat(I,1)=X(Iat)
163 1 pfleura2
     val_zmat(I,2)=Y(Iat)
164 1 pfleura2
     val_zmat(I,3)=Z(Iat)
165 1 pfleura2
     DejaFait(Iat)=.TRUE.
166 1 pfleura2
     idx_zmat(Iat)=I
167 1 pfleura2
  END DO
168 1 pfleura2
169 1 pfleura2
170 1 pfleura2
  izm=NCart
171 1 pfleura2
172 1 pfleura2
  ! We now calculate the connections
173 1 pfleura2
  Liaisons=0
174 1 pfleura2
  LiaisonsBis=0
175 1 pfleura2
176 1 pfleura2
  if (debug) THEN
177 1 pfleura2
     WRITE(*,*) "Liaisons initialized"
178 1 pfleura2
!     WRITE(*,*) 'Covalent radii used'
179 1 pfleura2
!     DO iat=1,na
180 1 pfleura2
!        i=atome(iat)
181 1 pfleura2
!        WRITE(*,*) Nom(I),Iat,r_cov(i),r_cov(i)*fact
182 1 pfleura2
!     END DO
183 1 pfleura2
  END IF
184 1 pfleura2
185 1 pfleura2
1003 FORMAT(1X,I4,' - ',25(I5))
186 1 pfleura2
187 1 pfleura2
  Call CalcCnct(na,atome,x,y,z,LIAISONS,r_cov,fact)
188 1 pfleura2
189 1 pfleura2
  if (debug) THEN
190 1 pfleura2
     WRITE(*,*) "Connections calculated"
191 1 pfleura2
     DO IL=1,na
192 1 pfleura2
        WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
193 1 pfleura2
     END DO
194 1 pfleura2
  END IF
195 1 pfleura2
196 1 pfleura2
  ! We get rid of bonds between cart atoms
197 1 pfleura2
  ! the trick is to do this on liaisons while keeping LiaisonsBis ok.
198 1 pfleura2
199 1 pfleura2
  LiaisonsIni=Liaisons
200 1 pfleura2
  LiaisonsBis=Liaisons
201 1 pfleura2
202 1 pfleura2
  DO I=1,NCart
203 1 pfleura2
     Iat=AtCart(I)
204 1 pfleura2
     if (debug) WRITE(*,'(20(1X,I3))') I,Iat,  &
205 1 pfleura2
          (LiaisonsIni(Iat,Itmp),ITmp=0,LiaisonsIni(Iat,0))
206 1 pfleura2
     DO J=1,LiaisonsIni(Iat,0)
207 1 pfleura2
        Jat=LiaisonsIni(Iat,J)
208 1 pfleura2
        IF (FCart(Jat)) THEN
209 1 pfleura2
           Call Annul(Liaisons,Iat,Jat)
210 1 pfleura2
           Call Annul(Liaisons,Jat,Iat)
211 1 pfleura2
           Call Annul(LiaisonsIni,Jat,Iat)
212 1 pfleura2
           Liaisons(Iat,0)=Liaisons(Iat,0)-1
213 1 pfleura2
           Liaisons(Jat,0)=Liaisons(Jat,0)-1
214 1 pfleura2
           LiaisonsIni(Jat,0)=LiaisonsIni(Jat,0)-1
215 1 pfleura2
        END IF
216 1 pfleura2
     END DO
217 1 pfleura2
  END DO
218 1 pfleura2
219 1 pfleura2
  if (debug) THEN
220 1 pfleura2
     WRITE(*,*) "Connections omitting bonds between cart atoms"
221 1 pfleura2
     DO IL=1,na
222 1 pfleura2
        WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
223 1 pfleura2
     END DO
224 1 pfleura2
  END IF
225 1 pfleura2
226 1 pfleura2
  ! We decompose the system into fragments, but we omit on purpuse
227 1 pfleura2
  ! fragments consisting only of frozen atoms
228 1 pfleura2
  ! Now, frozblock will contain the frozen atom indices in a given fragment
229 1 pfleura2
  Fragment=0
230 1 pfleura2
  NbAtFrag=0
231 1 pfleura2
  FrozBlock=0
232 1 pfleura2
233 1 pfleura2
  IdxFrag=0
234 1 pfleura2
  NbFrag=0
235 1 pfleura2
236 1 pfleura2
 Call Decomp_frag(na,liaisons,.NOT.Fcart,nbfrag,Fragment,NbAtFrag,FragAt)
237 1 pfleura2
238 1 pfleura2
 ! We compute FrozBlock now
239 1 pfleura2
  DO IFrag=1,NbFrag
240 1 pfleura2
     DO I=1,NbAtFrag(IFrag)
241 1 pfleura2
        Iat=FragAt(IFrag,I)
242 1 pfleura2
        IF (FCart(Iat)) THEN
243 1 pfleura2
           FrozBlock(IFrag,0)=FrozBlock(IFrag,0)+1
244 1 pfleura2
           FrozBlock(IFrag,FrozBlock(IFrag,0))=IAt
245 1 pfleura2
        END IF
246 1 pfleura2
     END DO
247 1 pfleura2
  END DO
248 1 pfleura2
249 1 pfleura2
  if (debug) THEN
250 1 pfleura2
     WRITE(*,*) 'I found ',NbFrag, 'fragments with respectively '
251 1 pfleura2
     WRITE(*,*) (NbAtFrag(I),I=1,NbFrag), 'atoms'
252 1 pfleura2
     WRITE(*,*) "Fragments atoms are now listed as F in the following"
253 1 pfleura2
     DO I=1,NbFrag
254 1 pfleura2
        WRITE(*,*) Na
255 1 pfleura2
        WRITE(*,*) 'Fragment ', i
256 1 pfleura2
        DO J=1,Na
257 1 pfleura2
           AtName=Nom(Atome(J))
258 1 pfleura2
           IF (Fragment(J).EQ.I) AtName='F'
259 1 pfleura2
           WRITE(*,'(1X,A5,3(1X,F10.6))') AtName,X(J),Y(J),Z(J)
260 1 pfleura2
        END DO
261 1 pfleura2
        !     WRITE(*,*) "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
262 1 pfleura2
     END DO
263 1 pfleura2
264 1 pfleura2
     DO I=1,NbFrag
265 1 pfleura2
        WRITE(*,*) 'Fragment ', i
266 1 pfleura2
        WRITE(*,'(A,20(1X,I3))') "FragAt:",(FragAt(I,j),j=1,NbAtFrag(I))
267 1 pfleura2
        WRITE(*,'(A,20(1X,I3))') "FrozBlock:",(FrozBlock(I,j),j=0,FrozBlock(I,0))
268 1 pfleura2
     END DO
269 1 pfleura2
  END IF
270 1 pfleura2
271 1 pfleura2
  NFroz=0
272 1 pfleura2
  DO IFrag=1,NbFrag
273 1 pfleura2
     If (FrozBlock(IFrag,0).NE.0) NFroz=NFroz+1
274 1 pfleura2
  END DO
275 1 pfleura2
276 1 pfleura2
  IF (DEBUG) WRITE(*,'(1X,A,I4,A,I4,A)') "Found ",NFroz," fragment(s) with cart atoms, out of",NbFrag," fragment(s)"
277 1 pfleura2
278 1 pfleura2
279 1 pfleura2
  if (debug) THEN
280 1 pfleura2
     WRITE(*,*) "Connections before adding fragments"
281 1 pfleura2
     DO IL=1,na
282 1 pfleura2
        IF (Liaisons(IL,0).NE.0)  WRITE(IOOUT,1003) Il,(LIAISONS(IL,JL),JL=0,NMaxL)
283 1 pfleura2
     END DO
284 1 pfleura2
  END IF
285 1 pfleura2
286 1 pfleura2
  ! We will now add the fragments containing non cart atoms.
287 1 pfleura2
  ! I am not sure that there is a best strategy to add them...
288 1 pfleura2
  ! so we add them in the order they were created :-(
289 1 pfleura2
  ! First only block with frozen atoms are added
290 1 pfleura2
  izm0=-1
291 1 pfleura2
  IFrag=1
292 1 pfleura2
  IZm=NCart
293 1 pfleura2
  FCaf=.FALSE.
294 1 pfleura2
  DO IFragFroz=1,NFroz
295 1 pfleura2
     DO WHILE (FrozBlock(IFrag,0).EQ.0)
296 1 pfleura2
        IFrag=IFrag+1
297 1 pfleura2
     END DO
298 1 pfleura2
     ! each atom has at least one frozen atom that will serve as the seed
299 1 pfleura2
     ! of this fragment.
300 1 pfleura2
     if (debug) WRITE(*,*) 'Adding fragment :', ifrag,'with',FrozBlock(IFrag,0),' frozen atoms'
301 1 pfleura2
     IF (FrozBlock(IFrag,0).EQ.1) THEN
302 1 pfleura2
        ! There is only one frozen atom, easy case ...
303 1 pfleura2
        Iat=FrozBlock(IFrag,1)
304 1 pfleura2
        if (debug) WRITE(*,*) 'Only frozen atom of frag',iat
305 1 pfleura2
        DejaFait(iat)=.TRUE.
306 1 pfleura2
        IdxCaFaire=2
307 1 pfleura2
        CaFaire(1)=iat
308 1 pfleura2
        CaFaire(2)=0
309 1 pfleura2
        IaFaire=1
310 1 pfleura2
        FCaf(Iat)=.TRUE.
311 1 pfleura2
        DO WHILE (CaFaire(IaFaire).NE.0)
312 1 pfleura2
           n1=CaFaire(IaFaire)
313 1 pfleura2
           I1=idx_zmat(n1)
314 1 pfleura2
           n2=ind_zmat(I1,2)
315 1 pfleura2
           n3=ind_zmat(I1,3)
316 1 pfleura2
317 1 pfleura2
        if (debug) WRITE(*,1003) n1,(LIAISONS(n1,JL),JL=0,NMaxL)
318 1 pfleura2
           DO I3=1,Liaisons(n1,0)
319 1 pfleura2
              IAt=Liaisons(n1,I3)
320 1 pfleura2
              ! PFL 2007.Apr.16 ->
321 1 pfleura2
              ! We add ALL atoms linked to n1 to CaFaire
322 1 pfleura2
              ! But we delete their link to n1
323 1 pfleura2
!! PFL 2007.Aug.28 ->
324 1 pfleura2
! re-add the test on FCaf in order not to put the same atom more than once in
325 1 pfleura2
! CaFaire array.
326 1 pfleura2
                 IF (.NOT.FCaf(Iat)) THEN
327 1 pfleura2
                    CaFaire(IdxCaFaire)=Iat
328 1 pfleura2
                    IdxCaFaire=IdxCaFaire+1
329 1 pfleura2
                    CaFaire(IdxCaFaire)=0
330 1 pfleura2
                    FCaf(Iat)=.TRUE.
331 1 pfleura2
                 END IF
332 1 pfleura2
!! <-- PFL 2007.Aug.28
333 1 pfleura2
              if (debug) WRITE(*,1003) Iat,(LIAISONS(Iat,JL),JL=0,NMaxL)
334 1 pfleura2
              Call Annul(Liaisons,Iat,n1)
335 1 pfleura2
              if (debug) WRITE(*,1003) Iat,(LIAISONS(Iat,JL),JL=0,NMaxL)
336 1 pfleura2
              Liaisons(iat,0)=Liaisons(Iat,0)-1
337 1 pfleura2
           if (debug) WRITE(*,*) 'Iafaire,n1,DejaFait,idxzmat,n2,n3',IaFaire,n1, &
338 1 pfleura2
                DejaFait(n1),idx_zmat(n1),n2,3
339 1 pfleura2
           if (debug) WRITE(*,*) 'Cafaire - 0', &
340 1 pfleura2
                (CaFaire(J),J=Iafaire,IdxCAfaire)
341 1 pfleura2
342 1 pfleura2
343 1 pfleura2
              ! <- PFL 2007.Apr.16
344 1 pfleura2
              IF (.NOT.DejaFait(Iat)) THEN
345 1 pfleura2
                 ! we add it to the Zmat
346 1 pfleura2
!                 WRITE(*,*) "coucou"
347 1 pfleura2
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
348 1 pfleura2
!                 WRITE(*,*) "coucou 2"
349 1 pfleura2
                 izm=izm+1
350 1 pfleura2
                 IF (izm.EQ.2) THEN
351 1 pfleura2
!                 WRITE(*,*) "coucou 3"
352 1 pfleura2
                    ind_zmat(izm,1)=IAt
353 1 pfleura2
                    ind_zmat(izm,2)=n1
354 1 pfleura2
                    ind_zmat(izm,3:5)=0
355 1 pfleura2
                    if (n2.EQ.-1) n2=Iat
356 1 pfleura2
                 END IF
357 1 pfleura2
358 1 pfleura2
!                 WRITE(*,*) "coucou 4"
359 1 pfleura2
                 if ((izm.GE.3).AND.(n2.EQ.-1)) THEN
360 1 pfleura2
!                 WRITE(*,*) "coucou 5"
361 1 pfleura2
!                 WRITE(*,*) "coucou 3 bis"
362 1 pfleura2
                       ! This atom is defined in cartesian coordinates
363 1 pfleura2
                       ! and has not been used to put a non cart atom.
364 1 pfleura2
                       ! we will find the closest atom linked to this one amongst the atoms
365 1 pfleura2
                       ! already stored in ind_zmat
366 1 pfleura2
                       d=1e9
367 1 pfleura2
                       JAt=-1
368 1 pfleura2
                       DO I=1,NCart
369 1 pfleura2
                          If (ind_zmat(I,1).NE.n1) THEN
370 1 pfleura2
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
371 1 pfleura2
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
372 1 pfleura2
                             ! we take only atoms that are not aligned
373 1 pfleura2
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
374 1 pfleura2
                                Jat=ind_zmat(I,1)
375 1 pfleura2
                                d=norm2
376 1 pfleura2
                             END IF
377 1 pfleura2
                          END IF
378 1 pfleura2
                       END DO
379 1 pfleura2
                       if (JAt.EQ.-1) THEN
380 1 pfleura2
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
381 1 pfleura2
                          STOP
382 1 pfleura2
                       END IF
383 1 pfleura2
                       n2=JAt
384 1 pfleura2
                 END IF! izm>=3 and n2.eq.-1
385 1 pfleura2
386 1 pfleura2
                 IF (izm.eq.3) THEN
387 1 pfleura2
                    ind_zmat(izm,1)=Iat
388 1 pfleura2
                    ind_zmat(izm,2)=n1
389 1 pfleura2
                    ind_zmat(izm,3)=n2
390 1 pfleura2
                 END IF  ! izm=3
391 1 pfleura2
392 1 pfleura2
                 IF (izm.GE.4) THEN
393 1 pfleura2
!                                     WRITE(*,*) "coucou 5 bis"
394 1 pfleura2
                    IF (n3.EQ.-1) THEN
395 1 pfleura2
                       ! This atom is defined in cartesian coordinates
396 1 pfleura2
                       ! and has not been used to put a non cart atom.
397 1 pfleura2
                       ! we will find the closest atom linked to this one amongst the atoms
398 1 pfleura2
                       ! already stored in ind_zmat
399 1 pfleura2
                       call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
400 1 pfleura2
                       d=1e9
401 1 pfleura2
                       JAt=-1
402 1 pfleura2
                       DO I=1,NCart
403 1 pfleura2
                          If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
404 1 pfleura2
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
405 1 pfleura2
                             SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
406 1 pfleura2
                             if (debug) WRITE(*,*) "1-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
407 1 pfleura2
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
408 1 pfleura2
                             ! we take only atoms that are not aligned
409 1 pfleura2
                             if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
410 1 pfleura2
                                Jat=ind_zmat(I,1)
411 1 pfleura2
                                d=norm3
412 1 pfleura2
                             END IF
413 1 pfleura2
                          END IF
414 1 pfleura2
                       END DO
415 1 pfleura2
                       if (JAt.EQ.-1) THEN
416 1 pfleura2
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
417 1 pfleura2
                          STOP
418 1 pfleura2
                       END IF
419 1 pfleura2
                       n3=JAt
420 1 pfleura2
                    END IF
421 1 pfleura2
                    !                    ind_zmat(I1,3)=n3
422 1 pfleura2
                    Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
423 1 pfleura2
                 END IF ! izm>=4
424 1 pfleura2
!                 WRITE(*,*) "coucou 6"
425 1 pfleura2
           if (debug) THEN
426 1 pfleura2
              WRITE(*,*) 'ind_zmat 0',izm
427 1 pfleura2
              DO Ip=1,NCart
428 1 pfleura2
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(ip,1:4)
429 1 pfleura2
                 END DO
430 1 pfleura2
              SELECT CASE (NCart)
431 1 pfleura2
                 CASE (1)
432 1 pfleura2
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
433 1 pfleura2
                    if (izm.GE.3) &
434 1 pfleura2
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
435 1 pfleura2
                    Ibeg=4
436 1 pfleura2
                 CASE (2)
437 1 pfleura2
                    WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
438 1 pfleura2
                    Ibeg=4
439 1 pfleura2
                 CASE DEFAULT
440 1 pfleura2
                    IBeg=NCart+1
441 1 pfleura2
                 END SELECT
442 1 pfleura2
              DO Ip=IBeg,izm
443 1 pfleura2
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
444 1 pfleura2
                      ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
445 1 pfleura2
              END DO
446 1 pfleura2
           END IF
447 1 pfleura2
448 1 pfleura2
                 idx_zmat(iat)=izm
449 1 pfleura2
!                 Call Annul(Liaisons,iat,n1)
450 1 pfleura2
!                 Liaisons(iat,0)=Liaisons(Iat,0)-1
451 1 pfleura2
                 !                  Call Annul(Liaisons,n1,iat)
452 1 pfleura2
                 n3=iat
453 1 pfleura2
                 DejaFait(Iat)=.TRUE.
454 1 pfleura2
              END IF
455 1 pfleura2
           END DO
456 1 pfleura2
           IaFaire=IaFaire+1
457 1 pfleura2
458 1 pfleura2
           if (debug) THEN
459 1 pfleura2
              WRITE(*,*) 'ind_zmat 5'
460 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
461 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
462 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
463 1 pfleura2
              DO Ip=4,izm
464 1 pfleura2
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
465 1 pfleura2
                      ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
466 1 pfleura2
              END DO
467 1 pfleura2
           END IF
468 1 pfleura2
469 1 pfleura2
        END Do              ! DO WHILE CaFaire
470 1 pfleura2
471 1 pfleura2
472 1 pfleura2
     ELSE     ! FrozBlock(I,0)==1
473 1 pfleura2
        if (debug) WRITE(*,*) 'IFrag=',IFrag,'Frozblock(I,0)>1',Frozblock(IFrag,0)
474 1 pfleura2
        ! We have many frozen atoms for one fragment. We will have to choose
475 1 pfleura2
        ! the first one, and then to decide when to swich...
476 1 pfleura2
        Iat=0
477 1 pfleura2
        I0=-1
478 1 pfleura2
        DO I=1,FrozBlock(IFrag,0)
479 1 pfleura2
           JAt=FrozBlock(IFrag,I)
480 1 pfleura2
           if (debug) WRITE(*,*) Jat, &
481 1 pfleura2
                (LiaisonsBis(Jat,Itmp),Itmp=0,LiaisonsBis(Jat,0))
482 1 pfleura2
           IF (LiaisonsBis(Jat,0).GT.I0) THEN
483 1 pfleura2
              I0=LiaisonsBis(Jat,0)
484 1 pfleura2
              Iat=Jat
485 1 pfleura2
           END IF
486 1 pfleura2
        END DO
487 1 pfleura2
        ! Iat is the starting frozen atom
488 1 pfleura2
        IF (debug) WRITE(*,*) 'Choosing atom ',iat,'as a starting atom for frag',ifrag
489 1 pfleura2
        DejaFait(iat)=.TRUE.
490 1 pfleura2
        IdxCaFaire=2
491 1 pfleura2
        CaFaire(1)=iat
492 1 pfleura2
        CaFaire(2)=0
493 1 pfleura2
        IaFaire=1
494 1 pfleura2
        FCaf(Iat)=.TRUE.
495 1 pfleura2
        DO WHILE (CaFaire(IaFaire).NE.0)
496 1 pfleura2
           n1=CaFaire(IaFaire)
497 1 pfleura2
           I1=idx_zmat(n1)
498 1 pfleura2
           n2=ind_zmat(I1,2)
499 1 pfleura2
           n3=ind_zmat(I1,3)
500 1 pfleura2
           if (debug) WRITE(*,*) 'Iafaire,n1,DejaFait,idxzmat,n2,n3',IaFaire,n1, &
501 1 pfleura2
                DejaFait(n1),idx_zmat(n1),n2,3
502 1 pfleura2
           if (debug) WRITE(*,*) 'Cafaire', &
503 1 pfleura2
                (CaFaire(J),J=Iafaire,IdxCAfaire)
504 1 pfleura2
505 1 pfleura2
506 1 pfleura2
           DO I3=1,Liaisons(n1,0)
507 1 pfleura2
              if (debug) WRITE(*,*) "Here, I3=",I3
508 1 pfleura2
              IAt=Liaisons(n1,I3)
509 1 pfleura2
              ! PFL 2007.Apr.16 ->
510 1 pfleura2
              ! We add ALL atoms linked to n1 to CaFaire
511 1 pfleura2
              ! But we delete their link to n1
512 1 pfleura2
!! PFL 2007.Aug.28 ->
513 1 pfleura2
! re-add the test on FCaf in order not to put the same atom more than once in
514 1 pfleura2
! CaFaire array.
515 1 pfleura2
                 IF (.NOT.FCaf(Iat)) THEN
516 1 pfleura2
                    CaFaire(IdxCaFaire)=Iat
517 1 pfleura2
                    IdxCaFaire=IdxCaFaire+1
518 1 pfleura2
                    CaFaire(IdxCaFaire)=0
519 1 pfleura2
                    FCaf(Iat)=.TRUE.
520 1 pfleura2
                 END IF
521 1 pfleura2
!! <-- PFL 2007.Aug.28
522 1 pfleura2
              Call Annul(Liaisons,Iat,n1)
523 1 pfleura2
              Liaisons(iat,0)=Liaisons(Iat,0)-1
524 1 pfleura2
              ! <- PFL 2007.Apr.16
525 1 pfleura2
              IF (.NOT.DejaFait(Iat)) THEN
526 1 pfleura2
                 ! we add it to the Zmat
527 1 pfleura2
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
528 1 pfleura2
                 izm=izm+1
529 1 pfleura2
                 IF (izm.EQ.2) THEN
530 1 pfleura2
                    ind_zmat(izm,1)=IAt
531 1 pfleura2
                    ind_zmat(izm,2)=n1
532 1 pfleura2
                    ind_zmat(izm,3:5)=0
533 1 pfleura2
                 ELSE
534 1 pfleura2
                    IF (n2.EQ.-1) THEN
535 1 pfleura2
                       ! This atom is defined in cartesian coordinates
536 1 pfleura2
                       ! and has not been used to put a non cart atom.
537 1 pfleura2
                       ! we will find the closest atom linked to this one amongst the atoms
538 1 pfleura2
                       ! already stored in ind_zmat
539 1 pfleura2
                       d=1e9
540 1 pfleura2
                       JAt=-1
541 1 pfleura2
                       DO I=1,NCart
542 1 pfleura2
                          If (ind_zmat(I,1).NE.n1) THEN
543 1 pfleura2
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
544 1 pfleura2
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
545 1 pfleura2
                             ! we take only atoms that are not aligned
546 1 pfleura2
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
547 1 pfleura2
                                Jat=ind_zmat(I,1)
548 1 pfleura2
                                d=norm2
549 1 pfleura2
                             END IF
550 1 pfleura2
                          END IF
551 1 pfleura2
                       END DO
552 1 pfleura2
                       if (JAt.EQ.-1) THEN
553 1 pfleura2
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
554 1 pfleura2
                          STOP
555 1 pfleura2
                       END IF
556 1 pfleura2
                       n2=JAt
557 1 pfleura2
                    END IF
558 1 pfleura2
                    !                    ind_zmat(I1,2)=n2
559 1 pfleura2
                    if (izm.EQ.3) THEN
560 1 pfleura2
                       ind_zmat(izm,1)=Iat
561 1 pfleura2
                       ind_zmat(izm,2)=n1
562 1 pfleura2
                       ind_zmat(izm,3)=n2
563 1 pfleura2
                    ELSE
564 1 pfleura2
                       IF (n3.EQ.-1) THEN
565 1 pfleura2
                          ! This atom is defined in cartesian coordinates
566 1 pfleura2
                          ! and has not been used to put a non cart atom.
567 1 pfleura2
                          ! we will find the closest atom linked to this one amongst the atoms
568 1 pfleura2
                          ! already stored in ind_zmat
569 1 pfleura2
                          call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
570 1 pfleura2
                          d=1e9
571 1 pfleura2
                          JAt=-1
572 1 pfleura2
                          DO I=1,NCart
573 1 pfleura2
                             If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
574 1 pfleura2
                                call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
575 1 pfleura2
                                SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
576 1 pfleura2
                             if (debug) WRITE(*,*) "2-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
577 1 pfleura2
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
578 1 pfleura2
579 1 pfleura2
                                ! we take only atoms that are not aligned
580 1 pfleura2
                                if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
581 1 pfleura2
                                   Jat=ind_zmat(I,1)
582 1 pfleura2
                                   d=norm3
583 1 pfleura2
                                END IF
584 1 pfleura2
                             END IF
585 1 pfleura2
                          END DO
586 1 pfleura2
                          if (JAt.EQ.-1) THEN
587 1 pfleura2
                             WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
588 1 pfleura2
                             STOP
589 1 pfleura2
                          END IF
590 1 pfleura2
                          n3=JAt
591 1 pfleura2
                       END IF
592 1 pfleura2
                       !                    ind_zmat(I1,3)=n3
593 1 pfleura2
                       Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
594 1 pfleura2
                    END IF
595 1 pfleura2
                 END IF
596 1 pfleura2
                 idx_zmat(iat)=izm
597 1 pfleura2
              !                  Call Annul(Liaisons,n1,iat)
598 1 pfleura2
                 n3=iat
599 1 pfleura2
                 DejaFait(Iat)=.TRUE.
600 1 pfleura2
                 if (debug) WRITE(*,*) "For no reason, I3=",I3
601 1 pfleura2
              END IF
602 1 pfleura2
           END DO
603 1 pfleura2
           IaFaire=IaFaire+1
604 1 pfleura2
605 1 pfleura2
           if (debug) THEN
606 1 pfleura2
              WRITE(*,*) 'ind_zmat 6',izm
607 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
608 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
609 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
610 1 pfleura2
              DO Ip=4,izm
611 1 pfleura2
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
612 1 pfleura2
                      ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
613 1 pfleura2
              END DO
614 1 pfleura2
           END IF
615 1 pfleura2
616 1 pfleura2
617 1 pfleura2
        END Do              ! DO WHILE CaFaire
618 1 pfleura2
619 1 pfleura2
620 1 pfleura2
     END IF  ! FrozBlock(I,0)==1
621 1 pfleura2
622 1 pfleura2
     IFrag=IFrag+1
623 1 pfleura2
624 1 pfleura2
  END DO                    ! Loop on all fragments
625 1 pfleura2
626 1 pfleura2
627 1 pfleura2
  DO IFrag=1,NbFrag
628 1 pfleura2
     IF (FrozBlock(IFrag,0).EQ.0) THEN
629 1 pfleura2
        if (debug) WRITE(*,*) 'Adding fragment :', ifrag,'with no frozen atoms'
630 1 pfleura2
        ! We have no frozen atoms for this fragment. We will have to choose
631 1 pfleura2
        ! the first atom to put.
632 1 pfleura2
        ! For now, we do not care of the 'central' atom (ie the one with
633 1 pfleura2
        ! no CP atoms...)
634 1 pfleura2
        ! We just take the atom that is the closest to the already placed
635 1 pfleura2
        ! atoms !
636 1 pfleura2
637 1 pfleura2
        I0=0
638 1 pfleura2
        I1=0
639 1 pfleura2
        d=1e9
640 1 pfleura2
        DO J=1,izm
641 1 pfleura2
           Jat=ind_zmat(J,1)
642 1 pfleura2
           DO I=1,NbAtfrag(IFrag)
643 1 pfleura2
              IAt=FragAt(IFrag,I)
644 1 pfleura2
              IF (.NOT.DejaFait(Iat)) THEN
645 1 pfleura2
                 Call vecteur(Iat,Jat,x,y,z,vx1,vy1,vz1,norm1)
646 1 pfleura2
                 IF (norm1.LT.d) THEN
647 1 pfleura2
                    I1=Jat
648 1 pfleura2
                    I0=Iat
649 1 pfleura2
                    d=Norm1
650 1 pfleura2
                 END IF
651 1 pfleura2
              END IF
652 1 pfleura2
           END DO
653 1 pfleura2
        END DO
654 1 pfleura2
655 1 pfleura2
        Iat=I0
656 1 pfleura2
        n1=I1
657 1 pfleura2
        I1=idx_zmat(n1)
658 1 pfleura2
        n2=ind_zmat(I1,2)
659 1 pfleura2
        n3=ind_zmat(I1,3)
660 1 pfleura2
661 1 pfleura2
        IF (debug) WRITE(*,*) Iat,' starting atom for frag',ifrag
662 1 pfleura2
663 1 pfleura2
664 1 pfleura2
        call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
665 1 pfleura2
        izm=izm+1
666 1 pfleura2
        IF (izm.EQ.2) THEN
667 1 pfleura2
           ind_zmat(izm,1)=IAt
668 1 pfleura2
           ind_zmat(izm,2)=n1
669 1 pfleura2
           ind_zmat(izm,3:5)=0
670 1 pfleura2
        ELSE
671 1 pfleura2
           IF (n2.EQ.-1) THEN
672 1 pfleura2
! This atom is defined in cartesian coordinates
673 1 pfleura2
! and has not been used to put a non cart atom.
674 1 pfleura2
! we will find the closest atom linked to this one amongst the atoms
675 1 pfleura2
! already stored in ind_zmat
676 1 pfleura2
              d=1e9
677 1 pfleura2
              JAt=-1
678 1 pfleura2
              DO I=1,NCart
679 1 pfleura2
                 If (ind_zmat(I,1).NE.n1) THEN
680 1 pfleura2
                    call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
681 1 pfleura2
                    SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
682 1 pfleura2
! we take only atoms that are not aligned
683 1 pfleura2
                    if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
684 1 pfleura2
                       Jat=ind_zmat(I,1)
685 1 pfleura2
                       d=norm2
686 1 pfleura2
                    END IF
687 1 pfleura2
                 END IF
688 1 pfleura2
              END DO
689 1 pfleura2
              if (JAt.EQ.-1) THEN
690 1 pfleura2
                 WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
691 1 pfleura2
                 STOP
692 1 pfleura2
              END IF
693 1 pfleura2
              n2=JAt
694 1 pfleura2
           END IF
695 1 pfleura2
!                    ind_zmat(I1,2)=n2
696 1 pfleura2
        END IF
697 1 pfleura2
        if (izm.EQ.3) THEN
698 1 pfleura2
           ind_zmat(izm,1)=Iat
699 1 pfleura2
           ind_zmat(izm,2)=n1
700 1 pfleura2
           ind_zmat(izm,3)=n2
701 1 pfleura2
        ELSE
702 1 pfleura2
           IF (n3.EQ.-1) THEN
703 1 pfleura2
! This atom is defined in cartesian coordinates
704 1 pfleura2
! and has not been used to put a non cart atom.
705 1 pfleura2
! we will find the closest atom linked to this one amongst the atoms
706 1 pfleura2
! already stored in ind_zmat
707 1 pfleura2
              call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
708 1 pfleura2
              d=1e9
709 1 pfleura2
              JAt=-1
710 1 pfleura2
              DO I=1,NCart
711 1 pfleura2
                 If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
712 1 pfleura2
                    call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
713 1 pfleura2
                    SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
714 1 pfleura2
                             if (debug) WRITE(*,*) "3-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
715 1 pfleura2
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
716 1 pfleura2
717 1 pfleura2
! we take only atoms that are not aligned
718 1 pfleura2
                    if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
719 1 pfleura2
                       Jat=ind_zmat(I,1)
720 1 pfleura2
                       d=norm3
721 1 pfleura2
                    END IF
722 1 pfleura2
                 END IF
723 1 pfleura2
              END DO
724 1 pfleura2
              if (JAt.EQ.-1) THEN
725 1 pfleura2
                 WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
726 1 pfleura2
                 STOP
727 1 pfleura2
              END IF
728 1 pfleura2
              n3=JAt
729 1 pfleura2
           END IF
730 1 pfleura2
           !                    ind_zmat(I1,3)=n3
731 1 pfleura2
           Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
732 1 pfleura2
        END IF
733 1 pfleura2
        idx_zmat(iat)=izm
734 1 pfleura2
735 1 pfleura2
736 1 pfleura2
        DejaFait(iat)=.TRUE.
737 1 pfleura2
        IdxCaFaire=2
738 1 pfleura2
        CaFaire(1)=iat
739 1 pfleura2
        CaFaire(2)=0
740 1 pfleura2
        IaFaire=1
741 1 pfleura2
        FCaf(Iat)=.TRUE.
742 1 pfleura2
        DO WHILE (CaFaire(IaFaire).NE.0)
743 1 pfleura2
           n1=CaFaire(IaFaire)
744 1 pfleura2
           if (debug) WRITE(*,*) 'Iafaire,n1',IaFaire,n1, &
745 1 pfleura2
                DejaFait(n1),idx_zmat(n1)
746 1 pfleura2
           if (debug) WRITE(*,*) 'Cafaire', &
747 1 pfleura2
                (CaFaire(J),J=Iafaire,IdxCAfaire)
748 1 pfleura2
           I1=idx_zmat(n1)
749 1 pfleura2
           n2=ind_zmat(I1,2)
750 1 pfleura2
           n3=ind_zmat(I1,3)
751 1 pfleura2
752 1 pfleura2
           DO I3=1,Liaisons(n1,0)
753 1 pfleura2
              IAt=Liaisons(n1,I3)
754 1 pfleura2
              ! PFL 2007.Apr.16 ->
755 1 pfleura2
              ! We add ALL atoms linked to n1 to CaFaire
756 1 pfleura2
              ! But we delete their link to n1
757 1 pfleura2
!! PFL 2007.Aug.28 ->
758 1 pfleura2
! re-add the test on FCaf in order not to put the same atom more than once in
759 1 pfleura2
! CaFaire array.
760 1 pfleura2
                 IF (.NOT.FCaf(Iat)) THEN
761 1 pfleura2
                    CaFaire(IdxCaFaire)=Iat
762 1 pfleura2
                    IdxCaFaire=IdxCaFaire+1
763 1 pfleura2
                    CaFaire(IdxCaFaire)=0
764 1 pfleura2
                    FCaf(Iat)=.TRUE.
765 1 pfleura2
                 END IF
766 1 pfleura2
!! <-- PFL 2007.Aug.28
767 1 pfleura2
              Call Annul(Liaisons,Iat,n1)
768 1 pfleura2
              Liaisons(iat,0)=Liaisons(Iat,0)-1
769 1 pfleura2
              ! <- PFL 2007.Apr.16
770 1 pfleura2
              IF (.NOT.DejaFait(Iat)) THEN
771 1 pfleura2
                 ! we add it to the Zmat
772 1 pfleura2
                 call vecteur(n1,iat,x,y,z,vx1,vy1,vz1,norm1)
773 1 pfleura2
                 izm=izm+1
774 1 pfleura2
                 IF (izm.EQ.2) THEN
775 1 pfleura2
                    ind_zmat(izm,1)=IAt
776 1 pfleura2
                    ind_zmat(izm,2)=n1
777 1 pfleura2
                    ind_zmat(izm,3:5)=0
778 1 pfleura2
                 ELSE
779 1 pfleura2
                    IF (n2.EQ.-1) THEN
780 1 pfleura2
! This atom is defined in cartesian coordinates
781 1 pfleura2
! and has not been used to put a non cart atom.
782 1 pfleura2
! we will find the closest atom linked to this one amongst the atoms
783 1 pfleura2
! already stored in ind_zmat
784 1 pfleura2
                       d=1e9
785 1 pfleura2
                       JAt=-1
786 1 pfleura2
                       DO I=1,NCart
787 1 pfleura2
                          If (ind_zmat(I,1).NE.n1) THEN
788 1 pfleura2
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx2,vy2,vz2,norm2)
789 1 pfleura2
                             SAng=Sinangle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
790 1 pfleura2
! we take only atoms that are not aligned
791 1 pfleura2
                             if ((norm2<=d).AND.(SAng.GT.0.09d0)) THEN
792 1 pfleura2
                                Jat=ind_zmat(I,1)
793 1 pfleura2
                                d=norm2
794 1 pfleura2
                             END IF
795 1 pfleura2
                          END IF
796 1 pfleura2
                       END DO
797 1 pfleura2
                       if (JAt.EQ.-1) THEN
798 1 pfleura2
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
799 1 pfleura2
                          STOP
800 1 pfleura2
                       END IF
801 1 pfleura2
                       n2=JAt
802 1 pfleura2
                    END IF
803 1 pfleura2
!                    ind_zmat(I1,2)=n2
804 1 pfleura2
                 END IF
805 1 pfleura2
                 if (izm.EQ.3) THEN
806 1 pfleura2
                    ind_zmat(izm,1)=Iat
807 1 pfleura2
                    ind_zmat(izm,2)=n1
808 1 pfleura2
                    ind_zmat(izm,3)=n2
809 1 pfleura2
                 ELSE
810 1 pfleura2
                    IF (n3.EQ.-1) THEN
811 1 pfleura2
! This atom is defined in cartesian coordinates
812 1 pfleura2
! and has not been used to put a non cart atom.
813 1 pfleura2
! we will find the closest atom linked to this one amongst the atoms
814 1 pfleura2
! already stored in ind_zmat
815 1 pfleura2
                       call vecteur(n1,n2,x,y,z,vx2,vy2,vz2,norm2)
816 1 pfleura2
                       d=1e9
817 1 pfleura2
                       JAt=-1
818 1 pfleura2
                       DO I=1,NCart
819 1 pfleura2
                          If ((ind_zmat(I,1).NE.n1).AND.(ind_zmat(I,1).NE.n2)) THEN
820 1 pfleura2
                             call vecteur(n1,ind_zmat(I,1),x,y,z,vx3,vy3,vz3,norm3)
821 1 pfleura2
                             SAng=Sinangle(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2)
822 1 pfleura2
                             if (debug) WRITE(*,*) "4-Trying izm,iat,n1,n2,n3,sang,d,norm3", &
823 1 pfleura2
                                  izm,iat,n1,n2,ind_zmat(I,1),sang,d,norm3
824 1 pfleura2
825 1 pfleura2
! we take only atoms that are not aligned
826 1 pfleura2
                             if ((norm3<=d).AND.(SAng.GT.0.09d0)) THEN
827 1 pfleura2
                                Jat=ind_zmat(I,1)
828 1 pfleura2
                                d=norm3
829 1 pfleura2
                             END IF
830 1 pfleura2
                          END IF
831 1 pfleura2
                       END DO
832 1 pfleura2
                       if (JAt.EQ.-1) THEN
833 1 pfleura2
                          WRITE(*,*) "All cart atoms are aligned. Case non treated for now. Contact PFL"
834 1 pfleura2
                          STOP
835 1 pfleura2
                       END IF
836 1 pfleura2
                       n3=JAt
837 1 pfleura2
                    END IF
838 1 pfleura2
!                    ind_zmat(I1,3)=n3
839 1 pfleura2
                    Call add2indzmat(na,izm,Iat,n1,n2,n3,ind_zmat,x,y,z)
840 1 pfleura2
                 END IF
841 1 pfleura2
                 idx_zmat(iat)=izm
842 1 pfleura2
                 n3=iat
843 1 pfleura2
                 DejaFait(Iat)=.TRUE.
844 1 pfleura2
              END IF
845 1 pfleura2
           END DO
846 1 pfleura2
           IaFaire=IaFaire+1
847 1 pfleura2
848 1 pfleura2
           if (debug) THEN
849 1 pfleura2
              WRITE(*,*) 'ind_zmat 7',izm
850 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(1,1)
851 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(2,1), ind_zmat(2,2)
852 1 pfleura2
              WRITE(*,'(1X,4(1X,I5))') ind_zmat(3,1), ind_zmat(3,2), ind_zmat(3,3)
853 1 pfleura2
              DO Ip=4,izm
854 1 pfleura2
                 WRITE(*,'(1X,4(1X,I5))') ind_zmat(Ip,1),  &
855 1 pfleura2
                      ind_zmat(Ip,2),ind_zmat(Ip,3), ind_zmat(Ip,4)
856 1 pfleura2
              END DO
857 1 pfleura2
           END IF
858 1 pfleura2
859 1 pfleura2
860 1 pfleura2
        END Do              ! DO WHILE CaFaire
861 1 pfleura2
     END IF                 ! FrozBlock(IFrag,0).EQ.0
862 1 pfleura2
863 1 pfleura2
  END DO                    ! Loop on all fragments without frozen atoms
864 1 pfleura2
865 1 pfleura2
  if (debug) THEN
866 1 pfleura2
     WRITE(*,*) Na+Izm
867 1 pfleura2
     WRITE(*,*) 'Done... :-(', izm
868 1 pfleura2
     DO J=1,Na
869 1 pfleura2
        WRITE(*,'(1X,A5,3(1X,F10.6))') Nom(Atome(J)),X(J),Y(J),Z(J)
870 1 pfleura2
     END DO
871 1 pfleura2
     DO I=1,Izm
872 1 pfleura2
        J=ind_zmat(I,1)
873 1 pfleura2
        WRITE(*,'(1X,A5,3(1X,F10.6))') 'X ',X(J),Y(J),Z(J)
874 1 pfleura2
     END DO
875 1 pfleura2
        IF (izm.NE.Na) THEN
876 1 pfleura2
           WRITE(*,*) "Atoms not done:"
877 1 pfleura2
           DO I=1,Na
878 1 pfleura2
              IF (.NOT.DejaFait(I)) WRITE(*,'(I5)',ADVANCE='NO') I
879 1 pfleura2
           END DO
880 1 pfleura2
        END IF
881 1 pfleura2
882 1 pfleura2
  END IF
883 1 pfleura2
884 1 pfleura2
885 1 pfleura2
  ! We have ind_zmat. We calculate val_zmat :-)
886 1 pfleura2
  if (debug) WRITE(*,*) "Calculating val_zmat"
887 1 pfleura2
888 1 pfleura2
  SELECT CASE (NCart)
889 1 pfleura2
     CASE (1)
890 1 pfleura2
        val_zmat(2,2)=0.d0
891 1 pfleura2
     val_zmat(2,3)=0.d0
892 1 pfleura2
     val_zmat(3,3)=0.d0
893 1 pfleura2
     n1=ind_zmat(2,1)
894 1 pfleura2
     n2=ind_zmat(2,2)
895 1 pfleura2
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
896 1 pfleura2
     val_zmat(2,1)=norm1
897 1 pfleura2
898 1 pfleura2
     n1=ind_zmat(3,1)
899 1 pfleura2
     n2=ind_zmat(3,2)
900 1 pfleura2
     n3=ind_zmat(3,3)
901 1 pfleura2
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
902 1 pfleura2
     val_zmat(3,1)=norm1
903 1 pfleura2
904 1 pfleura2
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
905 1 pfleura2
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
906 1 pfleura2
     val_zmat(3,2)=val
907 1 pfleura2
     IBeg=4
908 1 pfleura2
909 1 pfleura2
     CASE (2)
910 1 pfleura2
     val_zmat(3,3)=0.d0
911 1 pfleura2
912 1 pfleura2
     n1=ind_zmat(3,1)
913 1 pfleura2
     n2=ind_zmat(3,2)
914 1 pfleura2
     n3=ind_zmat(3,3)
915 1 pfleura2
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
916 1 pfleura2
     val_zmat(3,1)=norm1
917 1 pfleura2
918 1 pfleura2
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
919 1 pfleura2
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
920 1 pfleura2
     val_zmat(3,2)=val
921 1 pfleura2
     IBeg=4
922 1 pfleura2
  CASE DEFAULT
923 1 pfleura2
     Ibeg=NCart+1
924 1 pfleura2
  END SELECT
925 1 pfleura2
926 1 pfleura2
927 1 pfleura2
  DO i=IBeg,na
928 1 pfleura2
929 1 pfleura2
     n1=ind_zmat(i,1)
930 1 pfleura2
     n2=ind_zmat(i,2)
931 1 pfleura2
     n3=ind_zmat(i,3)
932 1 pfleura2
     n4=ind_zmat(i,4)
933 1 pfleura2
934 1 pfleura2
     if (debug) WRITE(*,*) "Doing i,n1,n2,n3,n4",i,n1,n2,n3,n4
935 1 pfleura2
     CALL vecteur(n2,n1,x,y,z,vx1,vy1,vz1,norm1)
936 1 pfleura2
937 1 pfleura2
     CALL vecteur(n2,n3,x,y,z,vx2,vy2,vz2,norm2)
938 1 pfleura2
     val=angle(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2)
939 1 pfleura2
940 1 pfleura2
     CALL vecteur(n3,n4,x,y,z,vx3,vy3,vz3,norm3)
941 1 pfleura2
     CALL produit_vect(vx1,vy1,vz1,norm1,vx2,vy2,vz2,norm2, &
942 1 pfleura2
          vx4,vy4,vz4,norm4)
943 1 pfleura2
     CALL produit_vect(vx3,vy3,vz3,norm3,vx2,vy2,vz2,norm2, &
944 1 pfleura2
          vx5,vy5,vz5,norm5)
945 1 pfleura2
946 1 pfleura2
     val_d=angle_d(vx4,vy4,vz4,norm4, vx5,vy5,vz5,norm5, &
947 1 pfleura2
          vx2,vy2,vz2,norm2)
948 1 pfleura2
949 1 pfleura2
     !     write(*,11) n1,n2,norm1,n3,val,n4,val_d
950 1 pfleura2
!11   format (2(1x,i3),1x,f8.4,2(1x,i3,1x,f8.3))
951 1 pfleura2
952 1 pfleura2
     val_zmat(i,1)=norm1
953 1 pfleura2
     val_zmat(i,2)=val
954 1 pfleura2
     val_zmat(i,3)=val_d
955 1 pfleura2
956 1 pfleura2
  END DO
957 1 pfleura2
958 1 pfleura2
  if (debug) THEN
959 1 pfleura2
     WRITE(*,*) 'DBG Cre_Zmat_Constr: ind_zmat'
960 1 pfleura2
     DO I=1,na
961 1 pfleura2
        WRITE(*,'(1X,5I5)') (ind_zmat(i,j),j=1,5)
962 1 pfleura2
     END DO
963 1 pfleura2
964 1 pfleura2
     WRITE(*,*) 'DBG Cre_Zmat_Constr: Full zmat'
965 1 pfleura2
     DO I=1,NCart
966 1 pfleura2
        WRITE(*,'(1X,I5,3(1X,F11.6))') ind_zmat(i,1),(val_zmat(i,j),j=1,3)
967 1 pfleura2
     END DO
968 1 pfleura2
     DO I=NCart+1,Na
969 1 pfleura2
        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)
970 1 pfleura2
     END DO
971 1 pfleura2
972 1 pfleura2
  END IF
973 1 pfleura2
974 1 pfleura2
  if (debug) WRITE(*,*) "Deallocate (FrozDist,Fragment, NbAtFrag,FragAt)"
975 1 pfleura2
  DEALLOCATE(FrozDist,Fragment, NbAtFrag,FragAt)
976 1 pfleura2
  if (debug) WRITE(*,*) "Deallocate FrozFrag, IdxFragAt, FrozBlock"
977 1 pfleura2
!  DEALLOCATE(FrozFrag, IdxFragAt, FrozBlock)
978 1 pfleura2
  DEALLOCATE(FrozFrag,FrozBlock)
979 1 pfleura2
  if (debug) WRITE(*,*) "Deallocate (DistFroz,Liaisons)"
980 1 pfleura2
  DEALLOCATE(DistFroz,Liaisons)
981 1 pfleura2
  if (debug) WRITE(*,*) "Deallocate (LiaisonsIni)"
982 1 pfleura2
  DEALLOCATE(LiaisonsIni)
983 1 pfleura2
  if (debug) WRITE(*,*) "Deallocate(CaFaire,DejaFait,FrozAt)"
984 1 pfleura2
  DEALLOCATE(CaFaire,DejaFait,FrozAt)
985 1 pfleura2
  if (debug) WRITE(*,*) "Deallocate LiaisonsBis"
986 1 pfleura2
  DEALLOCATE(LiaisonsBis)
987 1 pfleura2
  if (debug) WRITE(*,*) "Deallocate FCart,AtCart,FCaf"
988 1 pfleura2
  DEALLOCATE(FCart,AtCart,FCaf)
989 1 pfleura2
990 1 pfleura2
  if (debug) Call Header("Exiting Calc_mixed_frag")
991 1 pfleura2
992 1 pfleura2
END SUBROUTINE Calc_mixed_frag