root / src / Calc_mixed_frag.f90 @ 1
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 |