Statistiques
| Révision :

root / src / Calc_mixed_frag.f90 @ 2

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

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