root / utils / Xyz2Scan.f90 @ 6
Historique  Voir  Annoter  Télécharger (19,06 ko)
1 
program Xyz2irc 

2 
! This programs reads a XYZ file and converts it into distances, 
3 
! valence angle and dihedral angles. 
4 
! It prints them as a function of the irc distance... 
5 
! 
6 
! Input: name of the ROOT file of PAW 
7 
! it also needs a file call list which has the following structure: 
8 
! the first line gives the time when you want to start your analysis 
9 
! one line contains the type of the value you want to follow, it can be 
10 
! b for a Bond distance 
11 
! a for an angle 
12 
! d for a dihedral 
13 
! this descriptor is followed by the number of the atoms involved ! 
14 
! a typical file can be: 
15 
! 3. 
16 
! b 1 2 
17 
! b 2 3 
18 
! a 1 2 3 
19 
! 
20 
! Ouput: A files call Scan.dat 
21 
! wich contains in the first lines the input file (as a reminder) 
22 
! and then for each step the wanted values 
23 
! 
24 
! Second version also reads the energy (as to be written after E= on 
25 
! the comment line) 
26 
! 
27 
! Third version contains a new command: c for Center of Mass 
28 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
29 
! v 3.1 the c command now creates the center of mass... and allows 
30 
! people to do whatever they want with it... 
31 
! Syntax: c NbAt ListAt 
32 
!! 
33 
! v 3.2 
34 
! Added the p command that gives the oriented angle between two planes 
35 
! Syntax: p At1 At2 At3 At4 At5 At6 
36 
! At1, At2, At3 define the first plane 
37 
! At4, At5, At6 define the second plane 
38 
! How it works: gives the angles between the normal of the two planes, defined by 
39 
! the cross produc At2At1 x At2At3 etc. 
40 
!!!! 
41  
42  
43 
IMPLICIT NONE 
44 
INTEGER(4), PARAMETER :: maxnat=10000,MaxList=100 
45  
46 
CHARACTER(120) :: f1 
47 
REAL(8) :: geos(3,maxnat) 
48 
CHARACTER(33) :: fmt 
49 
CHARACTER(3) :: atoms(maxnat) 
50 
character(5) :: Type 
51 
CHARACTER(120) :: line 
52 
INTEGER(4) :: NbPrint 
53 
REAL(8) :: AU2PS,Pi 
54 
! Mass is not used for now. 
55 
! REAL(8) :: Mass(MaxNat) 
56 
REAL(8) :: Ener, Conv 
57 
INTEGER(4) :: At1,At2,At3,At4,At5,At6,IOOUT 
58 
INTEGER(4) :: IArg, I, NNN, Ng, J 
59  
60 
INTEGER(4) :: Nat,NbDist, NbAngle, NbDie,NbP,NbCOM 
61 
INTEGER(4) :: At1B(MaxList),At2B(MaxList) 
62 
INTEGER(4) :: At1A(MaxList),At2A(MaxList),At3A(MaxList) 
63 
INTEGER(4) :: At1D(MaxList),At2D(MaxList),At3D(MaxList),At4D(MaxList) 
64 
INTEGER(4) :: At1p(MaxList),At2p(MaxList),At3p(MaxList), & 
65 
At4p(MaxList),At5p(MaxList),At6p(MaxList) 
66 
INTEGER(4) :: AtCom(0:MaxList,MaxList) 
67 
REAL(8) :: VB(MaxList),VA(MaxList),VD(MaxList),VCOM(MaxList) 
68 
REAL(8) :: Vp(MaxList) 
69 
LOGICAL FExist 
70  
71 
COMMON /Indices/Nat,NbDist,NbAngle,NbDie,NbP,NbCom, At1B,At2B, & 
72 
At1A,At2A,At3A,At1D,At2D,At3D,At4D,AtCom, & 
73 
At1p,At2p,At3p,At4p,At5p,At6p 
74 
COMMON /Values/VB,VA,VD,Vp,VCom 
75 
COMMON /Const/AU2ps,Pi 
76  
77 
AU2PS=1./41341.37 
78 
NbPrint=100 
79 
Pi=dacos(1.d0) 
80 
IOOUT=13 
81 
Conv=1. 
82 
iarg=iargc() 
83 
if (iarg.lt.1) then 
84 
WRITE(*,*) "==============================================" 
85 
WRITE(*,*) "== ==" 
86 
WRITE(*,*) "== Xyz2Scan: Xyz file to scan file ==" 
87 
WRITE(*,*) "== ==" 
88 
WRITE(*,*) "==============================================" 
89 
WRITE(*,*) "Usage: xyz2scan XYZ_file " 
90 
WRITE(*,*) "The XYZ file should follow the 'usual' format:" 
91 
WRITE(*,*) "Number of atoms on the first line " 
92 
WRITE(*,*) "Comment on the second line" 
93 
WRITE(*,*) "The geometry is given in cartesian coordinates" 
94 
WRITE(*,*) "with atom symbol and three coordinates:" 
95 
WRITE(*,*) " C 0. 1. 0." 
96 
WRITE(*,*) "" 
97 
WRITE(*,*) " xyz2scan also needs a file called 'list' " 
98 
WRITE(*,*) "which has the following structure:" 
99 
WRITE(*,*) "Each line contains the type of the value you" 
100 
WRITE(*,*) "want to follow, it can be:" 
101 
WRITE(*,*) "  b for a Bond distance" 
102 
WRITE(*,*) "  a for an angle" 
103 
WRITE(*,*) "  d for a dihedral" 
104 
WRITE(*,*) "  p for the oriented angle between two planes" 
105 
WRITE(*,*) "This descriptor is followed by the number of" 
106 
WRITE(*,*) "the atoms involved" 
107 
WRITE(*,*) "One can also create 'barycenter' atoms that" 
108 
WRITE(*,*) "can then be used as normal atoms." 
109 
WRITE(*,*) "the descriptor is 'c' followed by the number" 
110 
WRITE(*,*) "of atoms and the list of atoms:" 
111 
WRITE(*,*) " c 2 1 5" 
112 
WRITE(*,*) "A typical file can be:" 
113 
WRITE(*,*) " b 1 2" 
114 
WRITE(*,*) " b 2 3" 
115 
WRITE(*,*) " a 1 2 3" 
116 
WRITE(*,*) " " 
117 
WRITE(*,*) "==============================================" 
118 
write(*,*) 'XYZ filename:' 
119 
read(*,'(a)') f1 
120 
else 
121 
call getarg(1,f1) 
122 
endif 
123  
124 
open(13,file='Scan.dat') 
125  
126 
INQUIRE(File='list',Exist=FExist) 
127 
if (.NOT.FExist) THEN 
128 
WRITE(*,*) "File *list* is missing" 
129 
STOP 
130 
END IF 
131  
132 
open(11,file=f1) 
133 
READ(11,*) nnn 
134 
close(11) 
135  
136 

137 
if (nnn.GT.MaxNat) THEN 
138 
WRITE(*,*) "Sorry but your system has too many atoms" 
139 
WRITE(*,*) "Change the value of MaxNat in the source file" 
140 
WRITE(*,*) "and then recompile." 
141 
WRITE(*,*) "For information, now MaxNat=",MaxNat 
142 
WRITE(*,*) "and you have ",nnn," atoms." 
143 
STOP 
144 
END IF 
145  
146 
Nat=nnn 
147  
148 
open(14,file='list') 
149 
Type="d" 
150 
DO WHILE (Type.ne."E") 
151 
CALL READLINE(14,Type,Line) 
152 
! WRITE(*,*) Line,Type 
153 
if (Type.eq."b") THEN 
154 
NbDist=NbDist+1 
155 
READ(Line,*) At1,At2 
156 
At1B(NbDist)=At1 
157 
At2B(NbDist)=At2 
158 
END IF 
159 
if (Type.eq."a") THEN 
160 
NbAngle=NbAngle+1 
161 
READ(Line,*) At1,At2,At3 
162 
At1A(NbAngle)=At1 
163 
At2A(NbAngle)=At2 
164 
At3A(NbAngle)=At3 
165 
END IF 
166 
if (Type.eq."d") THEN 
167 
NbDie=NbDie+1 
168 
READ(Line,*) At1,At2,At3,At4 
169 
At1D(NbDie)=At1 
170 
At2D(NbDie)=At2 
171 
At3D(NbDie)=At3 
172 
At4D(NbDie)=At4 
173 
! WRITE(13,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2, 
174 
! & Atoms(At3),At3,Atoms(At4),At4 
175 
! WRITE(*,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2, 
176 
! & Atoms(At3),At3,Atoms(At4),At4 
177  
178 
END IF 
179 
if (Type.eq."p") THEN 
180 
NbP=NbP+1 
181 
READ(Line,*) At1,At2,At3,At4,At5,At6 
182 
At1p(NbP)=At1 
183 
At2p(NbP)=At2 
184 
At3p(NbP)=At3 
185 
At4p(NbP)=At4 
186 
At5p(NbP)=At5 
187 
At6p(NbP)=At6 
188 
END IF 
189  
190 
if (Type.eq."c") THEN 
191 
NbCOM=NbCOM+1 
192 
READ(Line,*) At1,(AtCom(j,NbCOM),j=1,At1) 
193 
AtCom(0,NbCOM)=At1 
194 
Atoms(nat+NbCom)="G" 
195 
END IF 
196  
197 
END DO 
198  
199 
CLOSE(14) 
200  
201 
! We read one geoetry to get the atoms name 
202 
open(11,file=f1) 
203 
call rtraitem(11,nnn,ener,geos,atoms) 
204 
close(11) 
205  
206  
207 
! We write things 
208 
! First NbCom 
209 
if (NbCom.GE.1) THEN 
210 
WRITE(*,*) "# Added center of mass" 
211 
WRITE(IOOUT,*) "# Added center of mass" 
212 
DO J=1,NbCom 
213 
WRITE(IOOUT,'("# c ",I4,20(A3,I3))') AtCom(0,J), & 
214 
(Atoms(AtCom(i,J)),AtCom(i,J),i=1,At1) 
215 
WRITE(*,'("# c ",I4,20(A3,I3))') AtCom(0,J), & 
216 
(Atoms(AtCom(i,J)),AtCom(i,J),i=1,At1) 
217 
END DO 
218 
END IF 
219  
220 
! Distances 
221 
if (NbDist.GE.1) THEN 
222 
WRITE(*,*) "# Bonds" 
223 
WRITE(IOOUT,*) "# Bonds" 
224 
DO J=1,NbDist 
225 
At1= At1B(J) 
226 
At2=At2B(J) 
227 
WRITE(IOOUT,'("# b ",2(A3,I3))') Atoms(At1),At1, & 
228 
Atoms(At2),At2 
229 
WRITE(*,'("# b ",2(A3,I3))') Atoms(At1),At1,Atoms(At2),At2 
230 
END DO 
231 
END IF 
232  
233 
! Angles 
234 
if (NbAngle.GE.1) THEN 
235 
WRITE(*,*) "# Angles" 
236 
WRITE(IOOUT,*) "# Angles" 
237 
DO J=1,NbAngle 
238 
At1= At1A(J) 
239 
At2= At2A(J) 
240 
At3= At3A(J) 
241 
WRITE(IOOUT,'("# a ",3(A3,I3))') Atoms(At1),At1,Atoms(At2), & 
242 
At2, Atoms(At3),At3 
243 
WRITE(*,'("# a ",3(A3,I3))') Atoms(At1),At1,Atoms(At2), & 
244 
At2, Atoms(At3),At3 
245  
246 
END DO 
247 
END IF 
248  
249  
250 
! Dihedrals 
251 
if (NbDie.GE.1) THEN 
252 
WRITE(*,*) "# Dihedrals" 
253 
WRITE(IOOUT,*) "# Dihedrals" 
254 
DO J=1,NbDie 
255 
At1= At1D(J) 
256 
At2= At2D(J) 
257 
At3= At3D(J) 
258 
At4= At4D(J) 
259 
WRITE(IOOUT,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2) & 
260 
,At2,Atoms(At3),At3,Atoms(At4),At4 
261 
WRITE(*,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2, & 
262 
Atoms(At3),At3,Atoms(At4),At4 
263 
END DO 
264 
END IF 
265  
266 
! Planar angles 
267 
if (NbP.GE.1) THEN 
268 
WRITE(*,*) "# Planes angles" 
269 
WRITE(IOOUT,*) "# Planes angles" 
270 
DO J=1,NbP 
271 
At1= At1p(J) 
272 
At2= At2p(J) 
273 
At3= At3p(J) 
274 
At4= At4p(J) 
275 
At5= At5p(J) 
276 
At6= At6p(J) 
277 
WRITE(IOOUT,'("# p ",8(A3,I3))') Atoms(At1),At1,Atoms(At2) & 
278 
,At2,Atoms(At3),At3,Atoms(At4),At4 & 
279 
,Atoms(At5),At5,Atoms(At6),At6 
280 
WRITE(*,'("# p ",8(A3,I3))') Atoms(At1),At1,Atoms(At2),At2, & 
281 
Atoms(At3),At3,Atoms(At4),At4 & 
282 
,Atoms(At5),At5,Atoms(At6),At6 
283 
END DO 
284 
END IF 
285  
286  
287 
fmt='( (1X,F12.5),1X,F15.6)' 
288 
write(fmt(2:4),'(i3)') NbDist+NbAngle+NbDie+NbP 
289 
! write(*,*) nat3 
290 
! write(*,*) 'fmt:',fmt 
291  
292 
ng=1 
293  
294 
open(11,file=f1) 
295  
296 
10 call rtraitem(11,nnn,ener,geos,atoms) 
297 
! WRITE(*,*) nnn 
298  
299 
if (nnn.gt.0) then 
300 
! We convert coordinates in a0 into Angstroem 
301 
! write(*,*) "Analyse !" 
302 
Call Analyse(geos) 
303  
304 
write(IOOUT,fmt) (VB(j)/Conv,j=1,NbDist), & 
305 
(VA(j),j=1,NbAngle), & 
306 
(VD(j),j=1,NbDie),(Vp(j),j=1,NbP),ener 
307 
write(*,fmt) (VB(j)/Conv,j=1,NbDist), & 
308 
(VA(j),j=1,NbAngle), & 
309 
(VD(j),j=1,NbDie),(Vp(j),j=1,NbP),ener 
310  
311 
ng=ng+1 
312 
goto 10 
313 
endif 
314 
WRITE(*,*) 'Found ',ng1,' geometries' 
315 
close(11) 
316 
end 
317  
318 
! 
319 
subroutine rtraitem(ifil,nnn,E,tab,atoms) 
320 
! implicit REAL(8) :: (ah,oz) 
321 
IMPLICIT NONE 
322 
CHARACTER(120) :: line 
323 
CHARACTER(3) :: Atoms(*) 
324 
INTEGER(4) :: nnn, IFil, Idx, I, J 
325 
REAL(8) :: Tab(3,*), E 
326  
327  
328 
read(ifil,*,err=99,end=99) nnn 
329 
read(ifil,'(A)') line 
330 
idx=index(line,'E') 
331 
! WRITE(*,*) 'idx,line',idx,line 
332 
if (idx/=0) THEN 
333 
line=line(idx+2:) 
334 
idx=index(line,"=") 
335 
if (idx/=0) line=line(idx+1:) 
336 
idx=index(line,":") 
337 
if (idx/=0) line=line(idx+1:) 
338 
read(line,*) E 
339 
ELSE 
340 
E=0. 
341 
END IF 
342 

343 
! write(*,*) 'coucou',line 
344 
do i=1,nnn 
345 
read(ifil,'(A)',err=99,end=99) Line 
346 
do while (line(1:1)==' ') 
347 
line=line(2:) 
348 
end do 
349 
atoms(i)=line(1:3) 
350 
! write(*,*) 'coucou atoms',atoms(i) 
351 
do while (line(1:1).ne.' ') 
352 
line=line(2:) 
353 
end do 
354 
! WRITE(*,*) 'coucou2:',i,line 
355 
read(line,*) (tab(j,i),j=1,3) 
356 
end do 
357 
! WRITE(*,*) 'coucou:',nnn,tab(1,1) 
358 
return 
359 
99 nnn=0 
360 
! WRITE(*,*) 'Erreur lecture',ifil,nnn 
361 
return 
362 
end subroutine rtraitem 
363  
364 
! 
365  
366 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
367 
! READLINE 
368 
! This subroutine reads a line for a file, and converts 
369 
! the first field into a character variable, and the rest into 4 integers 
370 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
371  
372 
SUBROUTINE READLINE(IOIN,Type,Line) 
373 

374 
IMPLICIT NONE 
375 
CHARACTER(5) :: Type 
376 
CHARACTER(120) :: LINE 
377 
INTEGER(4) :: i,IOIN 
378 
READ(IOIN,'(A120)',ERR=999,END=999) LINE 
379 
! WRITE(*,*) Line 
380 
DO WHILE (LINE(1:1).eq.' ') 
381 
LINE=LINE(2:) 
382 
END DO 
383  
384 
i=1 
385 
DO WHILE (LINE(i:i).ne.' ') 
386 
i=i+1 
387 
END DO 
388 

389 
if (i.ge.6) THEN 
390 
WRITE(*,*) 'Pb with READLINE:',LINE 
391 
GOTO 999 
392 
END IF 
393 
Type=LINE(1:i1) 
394 
LINE=LINE(i:120) // " 0 0 0 0" 
395  
396 
RETURN 
397 
999 Type="E" 
398 
END SUBROUTINE READLINE 
399  
400  
401  
402 
SUBROUTINE Analyse(geos) 
403  
404 
IMPLICIT NONE 
405 
INTEGER(4),PARAMETER :: MaxNat=10000,MaxList=100 
406  
407 
REAL(8) :: geos(3,maxnat) 
408 
REAL(8) :: AU2PS,Pi 
409 
INTEGER(4) :: i,j,k 
410 
REAL(8) :: V1x,V1y,V1z,V2x,V2y,V2z,V3x,V3y,V3z 
411 
REAL(8) :: d1,d2,ca,sa 
412 
REAL(8) :: V4x,v4y,v4z,v5x,v5y,v5z 
413 
REAL(8) :: COG(3) 
414  
415 
INTEGER(4) :: Nat,NbDist, NbAngle, NbDie,NbP,NbCOM 
416 
INTEGER(4) :: At1B(MaxList),At2B(MaxList) 
417 
INTEGER(4) :: At1A(MaxList),At2A(MaxList),At3A(MaxList) 
418 
INTEGER(4) :: At1D(MaxList),At2D(MaxList),At3D(MaxList), & 
419 
At4D(MaxList) 
420 
INTEGER(4) :: At1p(MaxList),At2p(MaxList),At3p(MaxList), & 
421 
At4p(MaxList),At5p(MaxList),At6p(MaxList) 
422 
INTEGER(4) :: AtCom(0:MaxList,MaxList) 
423 
REAL(8) :: VB(MaxList),VA(MaxList),VD(MaxList),VCOM(MaxList) 
424 
REAL(8) :: Vp(MaxList) 
425  
426 
LOGICAL :: Debug=.FALSE. 
427  
428 
COMMON /Indices/Nat,NbDist,NbAngle,NbDie,NbP,NbCom, At1B,At2B, & 
429 
At1A,At2A,At3A,At1D,At2D,At3D,At4D,AtCom, & 
430 
At1p,At2p,At3p,At4p,At5p,At6p 
431 
COMMON /Values/VB,VA,VD,Vp,VCom 
432 
COMMON /Const/AU2ps,Pi 
433  
434 
if (Debug) THEN 
435 
DO I=1,Nat 
436 
WRITE(*,'(1X,I3,3(1X,F15.6))') i,(geos(j,i),j=1,3) 
437 
END DO 
438 
END IF 
439 
! First, we create the Centre of Mass atoms 
440 
DO i=1,NbCOM 
441 
COG(1)=0. 
442 
COG(2)=0. 
443 
COG(3)=0. 
444 
DO j=1,AtCOm(0,i) 
445 
DO k=1,3 
446 
COG(k)=COG(k)+geos(k,AtCom(j,i)) 
447 
END DO 
448 
DO k=1,3 
449 
COG(k)=COG(k)/AtCOM(0,i) 
450 
geos(k,Nat+i)=COG(k) 
451 
END DO 
452 
END DO 
453 
END DO 
454  
455 
DO i=1,NbDist 
456 
VB(i)=sqrt((geos(1,At1B(i))geos(1,At2B(i)))**2+ & 
457 
(geos(2,At1B(i))geos(2,At2B(i)))**2+ & 
458 
(geos(3,At1B(i))geos(3,At2B(i)))**2) 
459 
END DO 
460 
DO i=1,NbAngle 
461 
v1x=geos(1,At1A(i))geos(1,At2A(i)) 
462 
v1y=geos(2,At1A(i))geos(2,At2A(i)) 
463 
v1z=geos(3,At1A(i))geos(3,At2A(i)) 
464 
d1=sqrt(v1x**2+v1y**2+v1z**2) 
465 
v2x=geos(1,At3A(i))geos(1,At2A(i)) 
466 
v2y=geos(2,At3A(i))geos(2,At2A(i)) 
467 
v2z=geos(3,At3A(i))geos(3,At2A(i)) 
468 
d2=sqrt(v2x**2+v2y**2+v2z**2) 
469 
VA(i)=acos((v1x*v2x+v1y*v2y+v1z*v2z)/(d1*d2))*180./Pi 
470 
END DO 
471 
DO i=1,NbDie 
472 
! WRITE(*,*) At1D(i),At2D(i),At3D(i),At4D(i) 
473 
! WRITE(*,*) geos(1,At1D(i)),geos(2,At1D(i)),geos(3,At1D(i)) 
474 
! WRITE(*,*) geos(1,At2D(i)),geos(2,At2D(i)),geos(3,At2D(i)) 
475 
v1x=geos(1,At1D(i))geos(1,At2D(i)) 
476 
v1y=geos(2,At1D(i))geos(2,At2D(i)) 
477 
v1z=geos(3,At1D(i))geos(3,At2D(i)) 
478 
v2x=geos(1,At3D(i))geos(1,At2D(i)) 
479 
v2y=geos(2,At3D(i))geos(2,At2D(i)) 
480 
v2z=geos(3,At3D(i))geos(3,At2D(i)) 
481 
v3x=geos(1,At4D(i))geos(1,At3D(i)) 
482 
v3y=geos(2,At4D(i))geos(2,At3D(i)) 
483 
v3z=geos(3,At4D(i))geos(3,At3D(i)) 
484 

485 
v4x=v1y*v2zv1z*v2y 
486 
v4y=v1z*v2xv1x*v2z 
487 
v4z=v1x*v2yv1y*v2x 
488 
d1=sqrt(v4x**2+v4y**2+v4z**2) 
489 
v5x=v2y*v3z+v2z*v3y 
490 
v5y=v2z*v3x+v2x*v3z 
491 
v5z=v2x*v3y+v2y*v3x 
492 
d2=sqrt(v5x**2+v5y**2+v5z**2) 
493 
if (d1<=1e12) THEN 
494 
WRITE(*,*) "WARNING: Dihedral, d1=0" 
495 
ca=1. 
496 
sa=0. 
497 
END IF 
498 
if (d2<=1e12) THEN 
499 
WRITE(*,*) "WARNING: Dihedral, d2=0" 
500 
ca=1. 
501 
sa=0. 
502 
END IF 
503 
ca=(v4x*v5x+v4y*v5y+v4z*v5z)/(d1*d2) 
504 
sa=v1x*v5x+v1y*v5y+v1z*v5z 
505 
if (abs(ca)>1.) ca=sign(1.d0,ca) 
506 
VD(i)=acos(ca)*180./Pi 
507 
if (sa<0.) VD(i)=VD(i) 
508 
! WRITE(*,*) "Dihe",v5x,v5y,v5z,v4x,v4y,v4z,d1,d2, 
509 
! &(v4x*v5x+v4y*v5y+v4z*v5z)/(d1*d2),pi 
510 
! WRITE(*,*) "Dihe ca,sa,d1,d2",ca,sa,d1,d2,acos(ca) 
511  
512 

513 
!!!!!!!!! Another solution, more elegant ? 
514 
! norm2=sqrt(v2x**2+v2y**2+v2z**2) 
515 
! sa=(v4x*(v5y*v2zv5z*v2y) 
516 
! * v4y*(v5x*v2zv5z*v2x) 
517 
! * +v4z*(v5x*v2yv5y*v2x)) 
518 
! * /(d1*norm2*d2) 
519 
! angle_d=datan2(sa,ca)*180./Pi 
520 
! WRITE(*,*) sa,ca,angle_d,d1,d2,norm2 
521 
! WRITE(*,*) VD(i),angle_d 
522 
END DO 
523  
524 
DO i=1,NbP 
525 
! v1= At2At1 
526 
v1x=geos(1,At1p(i))geos(1,At2p(i)) 
527 
v1y=geos(2,At1p(i))geos(2,At2p(i)) 
528 
v1z=geos(3,At1p(i))geos(3,At2p(i)) 
529 
! v2=At2At3 
530 
v2x=geos(1,At3p(i))geos(1,At2p(i)) 
531 
v2y=geos(2,At3p(i))geos(2,At2p(i)) 
532 
v2z=geos(3,At3p(i))geos(3,At2p(i)) 
533 
! v4 = v1 x v2 
534 
v4x=v1y*v2zv1z*v2y 
535 
v4y=v1z*v2xv1x*v2z 
536 
v4z=v1x*v2yv1y*v2x 
537 
d1=sqrt(v4x**2+v4y**2+v4z**2) 
538  
539 
! v3= At5At4 
540 
v3x=geos(1,At4p(i))geos(1,At5p(i)) 
541 
v3y=geos(2,At4p(i))geos(2,At5p(i)) 
542 
v3z=geos(3,At4p(i))geos(3,At5p(i)) 
543 
! v2=At5At6 
544 
v2x=geos(1,At6p(i))geos(1,At5p(i)) 
545 
v2y=geos(2,At6p(i))geos(2,At5p(i)) 
546 
v2z=geos(3,At6p(i))geos(3,At5p(i)) 
547  
548 
! v5 = v3 x v2 
549 
v5x=v3y*v2zv3z*v2y 
550 
v5y=v3z*v2xv3x*v2z 
551 
v5z=v3x*v2yv3y*v2x 
552 
d2=sqrt(v5x**2+v5y**2+v5z**2) 
553  
554 
ca=(v4x*v5x+v4y*v5y+v4z*v5z)/(d1*d2) 
555 
sa=v1x*v5x+v1y*v5y+v1z*v5z 
556 
Vp(i)=acos(ca)*180./Pi 
557 
if (sa<0.) Vp(i)=Vp(i) 
558 
! WRITE(*,*) "Dihe",v5x,v5y,v5z,v4x,v4y,v4z,d1,d2, 
559 
! &(v4x*v5x+v4y*v5y+v4z*v5z)/(d1*d2),pi 
560 
!!!!!!!!! Another solution, more elegant ? 
561 
! See Dihedral routine 
562 
END DO 
563  
564  
565 
END SUBROUTINE Analyse 