root / utils / Xyz2Scan.f @ 3
History  View  Annotate  Download (18.8 kB)
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 maxnat,MaxList 
45 
Parameter (MaxNat=10000,MaxList=100) 
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), 
64 
$ At4D(MaxList) 
65 
INTEGER*4 At1p(MaxList),At2p(MaxList),At3p(MaxList), 
66 
$ At4p(MaxList),At5p(MaxList),At6p(MaxList) 
67 
INTEGER*4 AtCom(0:MaxList,MaxList) 
68 
REAL*8 VB(MaxList),VA(MaxList),VD(MaxList),VCOM(MaxList) 
69 
REAL*8 Vp(MaxList) 
70 
LOGICAL FExist 
71  
72 
c INTEGER*4 iargc 
73 
c external iargc 
74  
75 
COMMON /Indices/Nat,NbDist,NbAngle,NbDie,NbP,NbCom, At1B,At2B, 
76 
& At1A,At2A,At3A,At1D,At2D,At3D,At4D,AtCom, 
77 
& At1p,At2p,At3p,At4p,At5p,At6p 
78 
COMMON /Values/VB,VA,VD,Vp,VCom 
79 
COMMON /Const/AU2ps,Pi 
80  
81 
AU2PS=1./41341.37 
82 
NbPrint=100 
83 
Pi=dacos(1.d0) 
84 
IOOUT=13 
85 
Conv=1. 
86 
iarg=iargc() 
87 
if (iarg.lt.1) then 
88 
WRITE(*,*) "==============================================" 
89 
WRITE(*,*) "== ==" 
90 
WRITE(*,*) "== Xyz2Scan: Xyz file to scan file ==" 
91 
WRITE(*,*) "== ==" 
92 
WRITE(*,*) "==============================================" 
93 
WRITE(*,*) "Usage: xyz2scan XYZ_file " 
94 
WRITE(*,*) "The XYZ file should follow the 'usual' format:" 
95 
WRITE(*,*) "Number of atoms on the first line " 
96 
WRITE(*,*) "Comment on the second line" 
97 
WRITE(*,*) "The geometry is given in cartesian coordinates" 
98 
WRITE(*,*) "with atom symbol and three coordinates:" 
99 
WRITE(*,*) " C 0. 1. 0." 
100 
WRITE(*,*) "" 
101 
WRITE(*,*) " xyz2scan also needs a file called 'list' " 
102 
WRITE(*,*) "which has the following structure:" 
103 
WRITE(*,*) "Each line contains the type of the value you" 
104 
WRITE(*,*) "want to follow, it can be:" 
105 
WRITE(*,*) "  b for a Bond distance" 
106 
WRITE(*,*) "  a for an angle" 
107 
WRITE(*,*) "  d for a dihedral" 
108 
WRITE(*,*) "  p for the oriented angle between two planes" 
109 
WRITE(*,*) "This descriptor is followed by the number of" 
110 
WRITE(*,*) "the atoms involved" 
111 
WRITE(*,*) "One can also create 'barycenter' atoms that" 
112 
WRITE(*,*) "can then be used as normal atoms." 
113 
WRITE(*,*) "the descriptor is 'c' followed by the number" 
114 
WRITE(*,*) "of atoms and the list of atoms:" 
115 
WRITE(*,*) " c 2 1 5" 
116 
WRITE(*,*) "A typical file can be:" 
117 
WRITE(*,*) " b 1 2" 
118 
WRITE(*,*) " b 2 3" 
119 
WRITE(*,*) " a 1 2 3" 
120 
WRITE(*,*) " " 
121 
WRITE(*,*) "==============================================" 
122 
write(*,*) 'XYZ filename:' 
123 
read(*,'(a)') f1 
124 
else 
125 
call getarg(1,f1) 
126 
endif 
127  
128 
open(13,file='Scan.dat') 
129  
130 
INQUIRE(File='list',Exist=FExist) 
131 
if (.NOT.FExist) THEN 
132 
WRITE(*,*) "File *list* is missing" 
133 
STOP 
134 
END IF 
135  
136 
open(11,file=f1) 
137 
READ(11,*) nnn 
138 
close(11) 
139  
140 

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

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

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

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

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

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