root / utils / Xyz2Path.f90 @ 8
Historique  Voir  Annoter  Télécharger (29 ko)
1 
program Xyz2Path 

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 XYZ File 
7 
! it also needs a file call list which has the following structure: 
8 
! one line contains the type of the value you want to follow, it can be 
9 
! b for a Bond distance 
10 
! a for an angle 
11 
! d for a dihedral 
12 
! this descriptor is followed by the number of the atoms involved ! 
13 
! a typical file can be: 
14 
! b 1 2 
15 
! b 2 3 
16 
! a 1 2 3 
17 
! 
18 
! Ouput: A files call Scan.dat 
19 
! wich contains in the first lines the input file (as a reminder) 
20 
! and then for each step the wanted values 
21 
! 
22 
! Second version also reads the energy (as to be written after E= on 
23 
! the comment line) 
24 
! 
25 
! Third version contains a new command: c for Center of Mass 
26 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
27 
! v 3.1 the c command now creates the center of mass... and allows 
28 
! people to do whatever they want with it... 
29 
! Syntax: c NbAt ListAt 
30  
31  
32 
IMPLICIT NONE 
33  
34 
INCLUDE "Xyz2Path.param" 
35  
36 
character(40) :: f1 
37 
REAL(8) :: geos(3,maxnat), geos1(3,maxnat) 
38 
character(33) :: fmt 
39 
character(3) :: atoms(maxnat) 
40 
character(5) :: Type 
41 
Character(120) :: line 
42 
INTEGER(4) :: NbPrint 
43 
REAL(8) :: AU2PS,Pi 
44 
REAL(8) :: Mass(MaxNat), Ener, Conv, Ds, s 
45 
INTEGER(4) :: At1,At2,At3,At4,IOOUT,Iat 
46 
INTEGER(4) :: IArg, I, NNN, Ng, J 
47  
48 
INTEGER(4) :: Nat,NbDist, NbAngle, NbDie,NbCOM 
49 
INTEGER(4) :: At1B(MaxNat),At2B(MaxNat) 
50 
INTEGER(4) :: At1A(MaxNat),At2A(MaxNat),At3A(MaxNat) 
51 
INTEGER(4) :: At1D(MaxNat),At2D(MaxNat),At3D(MaxNat),At4D(MaxNat) 
52 
INTEGER(4) :: AtCom(0:MaxNat,MaxNat) 
53 
REAL(8) :: VB(MaxNat),VA(MaxNat),VD(MaxNat),VCOM(MaxNat) 
54  
55 
REAL(8) :: MRot(3,3), Rmsd 
56 
LOGICAL FExist,FRot,FAlign,Debug,DebugRMSD 
57  
58 
INTEGER(4) :: ConvertNumAt 
59 
external ConvertNumAt 
60  
61 
COMMON /Indices/Nat,NbDist,NbAngle,NbDie,NbCom, At1B,At2B, & 
62 
At1A,At2A,At3A,At1D,At2D,At3D,At4D,AtCom 
63 
COMMON /Values/VB,VA,VD,VCom 
64 
COMMON /Const/AU2ps,Pi 
65 

66  
67  
68 
REAL(8) :: MassAt(0:86)=(/0.0D0,1.0078D0, 4.0026D0, & 
69 
7.0160D0, 9.0122D0,11.0093D0, & 
70 
12.0000D0,14.0031D0,15.9949D0,18.9984D0,19.9924D0, & 
71 
22.9898D0,23.9850D0,26.9815D0, & 
72 
27.9769D0,30.9738D0,31.9721D0,34.9688D0,39.9624D0, & 
73 
39.0983D0,40.08D0, & 
74 
44.9559D0, 47.88D0, 50.9415D0, 51.996D0, 54.9380D0, & 
75 
55.847D0, 58.9332D0, 58.69D0, 63.546D0, 65.39D0, & 
76 
69.72D0,72.59D0,74.9216D0,78.96D0,79.904D0,83.80D0, & 
77 
85.4678D0,87.62D0,88.9059D0,91.224D0,92.9064D0, & 
78 
95.94D0,98D0,101.07D0,102.906D0,106.42D0,107.868D0,112.41D0, & 
79 
114.82D0,118.71D0,121.75D0,127.60D0,126.905D0,131.29D0, & 
80 
! 6 'CS','BA', 
81 
132.905D0,137.34D0, & 
82 
! 6 'LA', 
83 
! 'CE','PR','ND','PM','SM','EU','GD', 
84 
! 'TB','DY','HO', 'ER','TM','YB','LU', 
85 
138.91D0, & 
86 
140.12D0, 130.91D0, 144.24D0,147.D0,150.35D0, 151.96D0,157.25D0, & 
87 
158.924D0, 162.50D0, 164.93D0, 167.26D0,168.93D0,173.04D0,174.97D0, & 
88 
! 6 'HF','TA',' W','RE','OS','IR','PT', 
89 
! 'AU','HG', 
90 
! 6 'TL','PB','BI','PO','AT','RN'/ 
91 
178.49D0, 180.95D0, 183.85D0, 186.2D0, 190.2D0, 192.2D0, 195.09D0, & 
92 
196.97D0, 200.59D0, & 
93 
204.37D0, 207.19D0,208.98D0,210.D0,210.D0,222.D0 /) 
94  
95 
Debug=.FALSE. 
96 
DebugRMSD=.FALSE. 
97  
98 
AU2PS=1./41341.37 
99 
NbPrint=100 
100 
Pi=dacos(1.d0) 
101 
IOOUT=13 
102 
Conv=1. 
103 
FRot=.TRUE. 
104 
FAlign=.TRUE. 
105 
NbDist=0 
106 
NbAngle=0 
107 
NbDie=0 
108  
109 
iarg=command_argument_count() 
110 
if (iarg.lt.1) then 
111 
write(*,*) 'XYZ filename:' 
112 
read(*,'(a)') f1 
113 
else 
114 
call getarg(1,f1) 
115 
endif 
116  
117 
open(13,file='Scan.dat') 
118  
119 
INQUIRE(File='list',Exist=FExist) 
120 
if (.NOT.FExist) THEN 
121 
WRITE(*,*) "No file 'list', just printing Energy" 
122 
END IF 
123  
124 
open(11,file=f1) 
125 
call rtraitem(11,nnn,ener,geos1,atoms) 
126 
close(11) 
127  
128 
DO I=1,nnn 
129 
Iat=ConvertNumAt(atoms(I)) 
130 
Mass(I)=MassAt(Iat) 
131 
! write(*,*) I,Atoms(I),Mass(I),geos1(:,I) 
132 
END DO 
133  
134 
Nat=nnn 
135  
136 
if (FExist) THEN 
137 
open(14,file='list') 
138 
Type="d" 
139 
DO WHILE (Type.ne."E") 
140 
CALL READLINE(14,Type,Line) 
141 
! WRITE(*,*) Line,Type 
142 
if (Type.eq."b") THEN 
143 
NbDist=NbDist+1 
144 
READ(Line,*) At1,At2 
145 
At1B(NbDist)=At1 
146 
At2B(NbDist)=At2 
147 
WRITE(13,'("# b ",2(A3,I3))') Atoms(At1),At1,Atoms(At2),At2 
148 
WRITE(*,'("# b ",2(A3,I3))') Atoms(At1),At1,Atoms(At2),At2 
149 
END IF 
150 
if (Type.eq."a") THEN 
151 
NbAngle=NbAngle+1 
152 
READ(Line,*) At1,At2,At3 
153 
At1A(NbAngle)=At1 
154 
At2A(NbAngle)=At2 
155 
At3A(NbAngle)=At3 
156 
WRITE(13,'("# a ",3(A3,I3))') Atoms(At1),At1,Atoms(At2), & 
157 
At2, Atoms(At3),At3 
158 
WRITE(*,'("# a ",3(A3,I3))') Atoms(At1),At1,Atoms(At2), & 
159 
At2, Atoms(At3),At3 
160 
END IF 
161 
if (Type.eq."d") THEN 
162 
NbDie=NbDie+1 
163 
READ(Line,*) At1,At2,At3,At4 
164 
At1D(NbDie)=At1 
165 
At2D(NbDie)=At2 
166 
At3D(NbDie)=At3 
167 
At4D(NbDie)=At4 
168 
WRITE(13,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2, & 
169 
Atoms(At3),At3,Atoms(At4),At4 
170 
WRITE(*,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2, & 
171 
Atoms(At3),At3,Atoms(At4),At4 
172  
173 
END IF 
174 
if (Type.eq."c") THEN 
175 
NbCOM=NbCOM+1 
176 
READ(Line,*) At1,(AtCom(j,NbCOM),j=1,At1) 
177 
AtCom(0,NbCOM)=At1 
178 
Atoms(nat+NbCom)="G" 
179 
WRITE(13,'("# c ",I4,20(A3,I3))') At1, & 
180 
(Atoms(AtCom(i,NbCoM)),AtCom(i,NbCOM),i=1,At1) 
181 
WRITE(*,'("# c ",I4,20(A3,I3))') At1, & 
182 
(Atoms(AtCom(i,NbCoM)),AtCom(i,NbCOM),i=1,At1) 
183 
END IF 
184  
185 
END DO 
186 
END IF 
187  
188  
189  
190 
fmt='( (1X,F12.5),1X,F15.6)' 
191 
write(fmt(2:4),'(i3)') NbDist+NbAngle+NbDie+1 
192 
! write(*,*) nat3 
193 
! write(*,*) 'fmt:',fmt 
194  
195 
ng=1 
196  
197 
s=0. 
198  
199 
open(11,file=f1) 
200  
201 
10 call rtraitem(11,nnn,ener,geos,atoms) 
202 
! WRITE(*,*) nnn 
203 
if (nnn.gt.0) then 
204  
205 
if (debug) THEN 
206 
WRITE(*,*) "Debug ref geom: geos1" 
207 
DO I=1,nnn 
208 
WRITE(*,'(1X,I5,1X,F5.1,3(1X,F12.6))') I,mass(I),geos1(:,I) 
209 
END DO 
210 
END IF 
211  
212 
call CalcRmsd(nnn, geos1, geos,MRot,rmsd,FRot,FAlign,debugRMSD) 
213  
214 
if (debug) THEN 
215 
WRITE(*,*) "Debug Geos aligned on geos1,rmsd=",rmsd 
216 
DO I=1,nnn 
217 
WRITE(*,'(1X,I5,1X,F5.1,3(1X,F12.6))') I,mass(I),geos(:,I) 
218 
END DO 
219 
END IF 
220  
221 
ds=0. 
222 
DO I=1,nnn 
223 
DO J=1,3 
224 
ds=ds+mass(I)*(geos(J,I)geos1(J,I))**2 
225 
geos1(J,I)=Geos(J,I) 
226 
END DO 
227 
if (debug) write(*,*) I,ds 
228 
END DO 
229 
if (debug) WRITE(*,*) "Debug, ds",sqrt(ds) 
230 
s=s+sqrt(ds) 
231  
232 
! We convert coordinates in a0 into Angstroem 
233 
! write(*,*) "Analyse !" 
234 
if (FExist) THEN 
235 
Call Analyse(geos) 
236  
237 
write(IOOUT,fmt) s,(VB(j)/Conv,j=1,NbDist), & 
238 
(VA(j),j=1,NbAngle), & 
239 
(VD(j),j=1,NbDie),ener 
240 
write(*,fmt) s,(VB(j)/Conv,j=1,NbDist), & 
241 
(VA(j),j=1,NbAngle), & 
242 
(VD(j),j=1,NbDie),ener 
243 
ELSE 
244 
write(IOOUT,fmt) s,ener 
245 
write(*,fmt) s,ener 
246 
END IF 
247 
ng=ng+1 
248 
goto 10 
249 
endif 
250 
WRITE(*,*) 'Found ',ng1,' geometries' 
251 
close(11) 
252 
end 
253  
254 
! 
255 
subroutine rtraitem(ifil,nnn,E,tab,atoms) 
256 
! implicit REAL(8) :: (ah,oz) 
257 
IMPLICIT NONE 
258  
259 
character(120) :: line 
260 
character(3) :: Atoms(*) 
261 
INTEGER(4) :: nnn, IFil, Idx, I, J 
262 
REAL(8) :: Tab(3,*), E 
263  
264  
265 
read(ifil,*,err=99,end=99) nnn 
266 
read(ifil,'(A)') line 
267 
idx=index(line,'E') 
268 
! WRITE(*,*) 'idx,line',idx,line 
269 
if (idx/=0) THEN 
270 
line=line(idx+2:) 
271 
idx=index(line,"=") 
272 
if (idx/=0) line=line(idx+1:) 
273 
idx=index(line,":") 
274 
if (idx/=0) line=line(idx+1:) 
275 
read(line,*) E 
276 
ELSE 
277 
E=0. 
278 
END IF 
279 

280 
! write(*,*) 'coucou',line 
281 
do i=1,nnn 
282 
read(ifil,'(A)',err=99,end=99) Line 
283 
do while (line(1:1)==' ') 
284 
line=line(2:) 
285 
end do 
286 
atoms(i)=line(1:3) 
287 
! write(*,*) 'coucou atoms',atoms(i) 
288 
do while (line(1:1).ne.' ') 
289 
line=line(2:) 
290 
end do 
291 
! WRITE(*,*) 'coucou2:',i,line 
292 
read(line,*) (tab(j,i),j=1,3) 
293 
end do 
294 
! WRITE(*,*) 'coucou:',nnn,tab(1,1) 
295 
return 
296 
99 nnn=0 
297 
! WRITE(*,*) 'Erreur lecture',ifil,nnn 
298 
return 
299 
end 
300  
301 
! 
302  
303 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
304 
! READLINE 
305 
! This subroutine reads a line for a file, and converts 
306 
! the first field into a character variable, and the rest into 4 integers 
307 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
308  
309 
SUBROUTINE READLINE(IOIN,Type,Line) 
310 

311 
IMPLICIT NONE 
312 
CHARACTER(5) :: Type 
313 
CHARACTER(120) :: LINE 
314 
INTEGER(4) :: i,IOIN 
315 
READ(IOIN,'(A120)',ERR=999,END=999) LINE 
316 
! WRITE(*,*) Line 
317 
DO WHILE (LINE(1:1).eq.' ') 
318 
LINE=LINE(2:) 
319 
END DO 
320  
321 
i=1 
322 
DO WHILE (LINE(i:i).ne.' ') 
323 
i=i+1 
324 
END DO 
325 

326 
if (i.ge.6) THEN 
327 
WRITE(*,*) 'Pb with READLINE:',LINE 
328 
GOTO 999 
329 
END IF 
330 
Type=LINE(1:i1) 
331 
LINE=LINE(i:120) // " 0 0 0 0" 
332  
333 
RETURN 
334 
999 Type="E" 
335 
END 
336  
337  
338  
339 
SUBROUTINE Analyse(geos) 
340  
341 
IMPLICIT NONE 
342 
INCLUDE "Xyz2Path.param" 
343  
344 
REAL(8) :: geos(3,maxnat) 
345 
REAL(8) :: AU2PS,Pi 
346 
INTEGER(4) :: i,j,k 
347 
REAL(8) :: V1x,V1y,V1z,V2x,V2y,V2z,V3x,V3y,V3z 
348 
REAL(8) :: d1,d2,ca,sa 
349 
REAL(8) :: V4x,v4y,v4z,v5x,v5y,v5z 
350 
REAL(8) :: COG(3) 
351 
LOGICAL :: Debug=.FALSE. 
352 
INTEGER(4) :: Nat,NbDist, NbAngle, NbDie,NbCOM 
353 
INTEGER(4) :: At1B(MaxNat),At2B(MaxNat) 
354 
INTEGER(4) :: At1A(MaxNat),At2A(MaxNat),At3A(MaxNat) 
355 
INTEGER(4) :: At1D(MaxNat),At2D(MaxNat),At3D(MaxNat),At4D(MaxNat) 
356 
INTEGER(4) :: AtCom(0:MaxNat,MaxNat) 
357 
REAL(8) :: VB(MaxNat),VA(MaxNat),VD(MaxNat),VCOM(MaxNat) 
358  
359 
COMMON /Indices/Nat,NbDist,NbAngle,NbDie,NbCom, & 
360 
At1B,At2B,At1A,At2A,At3A,At1D,At2D,At3D,At4D,AtCom 
361 
COMMON /Values/VB,VA,VD,VCom 
362 
COMMON /Const/AU2ps,Pi 
363  
364 
if (debug) THEN 
365 
DO I=1,Nat 
366 
WRITE(*,'(1X,I3,3(1X,F15.6))') i,(geos(j,i),j=1,3) 
367 
END DO 
368 
END IF 
369  
370 
! First, we create the Centre of Mass atoms 
371 
DO i=1,NbCOM 
372 
COG(1)=0. 
373 
COG(2)=0. 
374 
COG(3)=0. 
375 
DO j=1,AtCOm(0,i) 
376 
DO k=1,3 
377 
COG(k)=COG(k)+geos(k,AtCom(j,i)) 
378 
END DO 
379 
END DO 
380 
DO k=1,3 
381 
COG(k)=COG(k)/AtCOM(0,i) 
382 
geos(k,Nat+i)=COG(k) 
383 
END DO 
384 
END DO 
385  
386  
387 
DO i=1,NbDist 
388 
VB(i)=sqrt((geos(1,At1B(i))geos(1,At2B(i)))**2+ & 
389 
(geos(2,At1B(i))geos(2,At2B(i)))**2+ & 
390 
(geos(3,At1B(i))geos(3,At2B(i)))**2) 
391 
END DO 
392 
DO i=1,NbAngle 
393 
v1x=geos(1,At1A(i))geos(1,At2A(i)) 
394 
v1y=geos(2,At1A(i))geos(2,At2A(i)) 
395 
v1z=geos(3,At1A(i))geos(3,At2A(i)) 
396 
d1=sqrt(v1x**2+v1y**2+v1z**2) 
397 
v2x=geos(1,At3A(i))geos(1,At2A(i)) 
398 
v2y=geos(2,At3A(i))geos(2,At2A(i)) 
399 
v2z=geos(3,At3A(i))geos(3,At2A(i)) 
400 
d2=sqrt(v2x**2+v2y**2+v2z**2) 
401 
VA(i)=acos((v1x*v2x+v1y*v2y+v1z*v2z)/(d1*d2))*180./Pi 
402 
END DO 
403 
DO i=1,NbDie 
404 
! WRITE(*,*) At1D(i),At2D(i),At3D(i),At4D(i) 
405 
! WRITE(*,*) geos(1,At1D(i)),geos(2,At1D(i)),geos(3,At1D(i)) 
406 
! WRITE(*,*) geos(1,At2D(i)),geos(2,At2D(i)),geos(3,At2D(i)) 
407 
v1x=geos(1,At1D(i))geos(1,At2D(i)) 
408 
v1y=geos(2,At1D(i))geos(2,At2D(i)) 
409 
v1z=geos(3,At1D(i))geos(3,At2D(i)) 
410 
v2x=geos(1,At3D(i))geos(1,At2D(i)) 
411 
v2y=geos(2,At3D(i))geos(2,At2D(i)) 
412 
v2z=geos(3,At3D(i))geos(3,At2D(i)) 
413 
v3x=geos(1,At4D(i))geos(1,At3D(i)) 
414 
v3y=geos(2,At4D(i))geos(2,At3D(i)) 
415 
v3z=geos(3,At4D(i))geos(3,At3D(i)) 
416 

417 
v4x=v1y*v2zv1z*v2y 
418 
v4y=v1z*v2xv1x*v2z 
419 
v4z=v1x*v2yv1y*v2x 
420 
d1=sqrt(v4x**2+v4y**2+v4z**2) 
421 
v5x=v2y*v3z+v2z*v3y 
422 
v5y=v2z*v3x+v2x*v3z 
423 
v5z=v2x*v3y+v2y*v3x 
424 
d2=sqrt(v5x**2+v5y**2+v5z**2) 
425 
ca=(v4x*v5x+v4y*v5y+v4z*v5z)/(d1*d2) 
426 
sa=v1x*v5x+v1y*v5y+v1z*v5z 
427 
VD(i)=acos(ca)*180./Pi 
428 
if (sa<0.) VD(i)=VD(i) 
429 
! WRITE(*,*) "Dihe",v5x,v5y,v5z,v4x,v4y,v4z,d1,d2, 
430 
! &(v4x*v5x+v4y*v5y+v4z*v5z)/(d1*d2),pi 
431 
!!!!!!!!! Another solution, more elegant ? 
432 
! norm2=sqrt(v2x**2+v2y**2+v2z**2) 
433 
! sa=(v4x*(v5y*v2zv5z*v2y) 
434 
! * v4y*(v5x*v2zv5z*v2x) 
435 
! * +v4z*(v5x*v2yv5y*v2x)) 
436 
! * /(d1*norm2*d2) 
437 
! angle_d=datan2(sa,ca)*180./Pi 
438 
! WRITE(*,*) sa,ca,angle_d,d1,d2,norm2 
439 
! WRITE(*,*) VD(i),angle_d 
440 
END DO 
441  
442 
END SUBROUTINE Analyse 
443  
444 
!C================================================================ 
445 
!C Convertit un nom d'atome (2 lettres) en nombre de masse (entier) 
446 
!C cette fonction a ete modifiee pour pouvoir convertir un nom 
447 
!C d'atome avec un numero a la fin... 
448 
!C================================================================ 
449  
450 
FUNCTION ConvertNumAt(ATOM) 
451  
452  
453 
IMPLICIT NONE 
454  
455 
INTEGER(4) :: I,Long,ConvertNumAt,IC 
456 
character(2) :: ATOM,L_Atom 
457 
character(3) :: ATOME 
458 
INTEGER(4), PARAMETER :: Max_Z=86 
459  
460  
461 
CHARACTER(2) :: Nom(0:max_Z)=(/ ' X',' H', 'HE', & 
462 
'LI','BE', ' B',' C',' N',' O',' F','NE', & 
463 
'NA','MG', 'AL','SI',' P',' S','CL','AR', & 
464 
' K','CA', & 
465 
'SC','TI',' V','CR','MN','FE','CO','NI','CU','ZN', & 
466 
'GA','GE','AS','SE','BR','KR', & 
467 
'RB','SR', & 
468 
' Y','ZR','NB','MO','TC','RU','RH','PD','AG','CD', & 
469 
'IN','SN','SB','TE',' I','XE', & 
470 
'CS','BA', & 
471 
'LA', & 
472 
'CE','PR','ND','PM','SM','EU','GD','TB','DY','HO', & 
473 
'ER','TM','YB','LU', & 
474 
'HF','TA',' W','RE','OS','IR','PT','AU','HG', & 
475 
'TL','PB','BI','PO','AT','RN'/) 
476  
477  
478 
!C Verifie qu'il n'y a que des lettres et des espaces dans ATOM 
479 
! WRITE(*,*) 'DBG CNVNUMAT, ATOM:',ATOM 
480  
481 
L_atom=Atom 
482 
IF (ATOM(1:1).LT.'A') L_ATOM(1:1)=' ' 
483 
IC=Ichar(ATOM(1:1)) 
484 
IF ((ic.le.123).AND.(ic.ge.97)) L_ATOM(1:1)=CHAr(IC32) 
485 
IF (ATOM(2:2).LT.'A') L_ATOM(2:2)=' ' 
486 
IC=Ichar(ATOM(2:2)) 
487 
IF ((ic.le.123).AND.(ic.ge.97)) L_ATOM(2:2)=CHAr(IC32) 
488 
!C Justifie le nom sur la droite (et non sur la gauche comme souvent...) 
489 
Long=INDEX(L_ATOM,' ')1 
490 
ATOME=' ' // L_ATOM 
491 
IF (Long.EQ.1) L_ATOM=ATOME(1:2) 
492 
! WRITE(*,*) 'DBG CNVNUMAT, L_ATOM:',L_ATOM 
493 
I=max_Z 
494 
DO WHILE ((nom(I).NE.L_ATOM) .AND. (I.GT.0)) 
495 
I=I1 
496 
END DO 
497 
ConvertNumAT=I 
498 
END 
499 
! This subroutine calculates RMSD using quaternions. 
500 
! It is based on the F90 routine bu E. Coutsias 
501 
! http://www.math.unm.edu/~vageli/homepage.html 
502 
! I (PFL) have just translated it, and I have changed the diagonalization 
503 
! subroutine. 
504 
! I also made some changes to make it suitable for Cart package. 
505 
! 
506 
! 
507 
! Copyright (C) 2004, 2005 Chaok Seok, Evangelos Coutsias and Ken Dill 
508 
! UCSF, Univeristy of New Mexico, Seoul National University 
509 
! Witten by Chaok Seok and Evangelos Coutsias 2004. 
510 

511 
! This library is free software; you can redistribute it and/or 
512 
! modify it under the terms of the GNU Lesser General Public 
513 
! License as published by the Free Software Foundation; either 
514 
! version 2.1 of the License, or (at your option) any later version. 
515 
! 
516 

517 
! This library is distributed in the hope that it will be useful, 
518 
! but WITHOUT ANY WARRANTY; without even the implied warranty of 
519 
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
520 
! Lesser General Public License for more details. 
521 
! 
522 

523 
! You should have received a copy of the GNU Lesser General Public 
524 
! License along with this library; if not, write to the Free Software 
525 
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA 
526 
! 
527  
528 
subroutine CalcRmsd(na, geom,geom2,U,rmsd,FRot,FAlign,debug) 
529 
! 
530 
! This subroutine calculates the least square rmsd of two coordinate 
531 
! sets coord1(3,n) and coord2(3,n) using a method based on quaternion. 
532 
! It then calculate the rotation matrix U and the centers of coord, and uses 
533 
! them to align the two molecules. 
534 
! 
535  
536  
537 
IMPLICIT NONE 
538  
539 
INCLUDE "Xyz2Path.param" 
540  
541 
INTEGER(4) :: IOIN, IOOUT, IOSCRN, IOPRNT, IOTAMP 
542 
COMMON/IODEFS/IOIN,IOOUT,IOSCRN,IOPRNT,IOTAMP 
543  
544 
INTEGER(4) :: na 
545 
REAL(8) :: geom(3,MaxNAt), Geom2(3,MaxNAt) 
546 
REAL(8) :: U(3,3), rmsd 
547 
LOGICAL FRot,FAlign,Debug 
548  
549 
REAL(8) :: Coord1(3,MaxNAt), Coord2(3,MaxNAt) 
550 
REAL(8) :: xc1,yc1,zc1, xc2,yc2,zc2 
551 

552  
553 
INTEGER(4) :: i, j, ia 
554 
REAL(8) :: x_norm, y_norm, lambda 
555 
REAL(8) :: Rmatrix(3,3) 
556 
REAL(8) :: S(4,4) 
557 
REAL(8) :: EigVec(4,4), EigVal(4) 
558  
559  
560  
561 
! calculate the barycenters, centroidal coordinates, and the norms 
562 
x_norm = 0.0d0 
563 
y_norm = 0.0d0 
564 
xc1=0. 
565 
yc1=0. 
566 
zc1=0. 
567 
xc2=0. 
568 
yc2=0. 
569 
zc2=0. 
570 
do ia=1,na 
571 
xc1=xc1+geom(1,ia) 
572 
xc2=xc2+geom2(1,ia) 
573 
yc1=yc1+geom(2,ia) 
574 
yc2=yc2+geom2(2,ia) 
575 
zc1=zc1+geom(3,ia) 
576 
zc2=zc2+geom2(3,ia) 
577 
! if (debug) WRITE(*,'(A,I5,4(1X,F10.4))') 'ia...',ia,x(ia), 
578 
! & x2(ia),xc1,xc2 
579 
END DO 
580 
xc1=xc1/dble(na) 
581 
yc1=yc1/dble(na) 
582 
zc1=zc1/dble(na) 
583 
xc2=xc2/dble(na) 
584 
yc2=yc2/dble(na) 
585 
zc2=zc2/dble(na) 
586  
587 
IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center1',xc1,yc1,zc1 
588 
IF (debug) WRITE(*,'(1X,A,3(1X,F10.4))') 'Center2',xc2,yc2,zc2 
589 
do i=1,na 
590 
Coord1(1,i)=geom(1,i)xc1 
591 
Coord1(2,i)=geom(2,i)yc1 
592 
Coord1(3,i)=geom(3,i)zc1 
593 
Coord2(1,i)=geom2(1,i)xc2 
594 
Coord2(2,i)=geom2(2,i)yc2 
595 
Coord2(3,i)=geom2(3,i)zc2 
596 
x_norm=x_norm+Coord1(1,i)**2+Coord1(2,i)**2+Coord1(3,i)**2 
597 
y_norm=y_norm+Coord2(1,i)**2+Coord2(2,i)**2+Coord2(3,i)**2 
598 
end do 
599  
600 
IF (debug) THEN 
601 
WRITE(*,*) "R matrix" 
602 
DO I=1,3 
603 
WRITE(*,*) (RMatrix(I,j),j=1,3) 
604 
END DO 
605 
END IF 
606  
607 
! calculate the R matrix 
608 
do i = 1, 3 
609 
do j = 1, 3 
610 
Rmatrix(i,j)=0. 
611 
do ia=1,na 
612 
Rmatrix(i,j) = Rmatrix(i,j)+Coord1(i,ia)*Coord2(j,ia) 
613 
END DO 
614 
end do 
615 
end do 
616  
617 
IF (debug) THEN 
618 
WRITE(*,*) "R matrix" 
619 
DO I=1,3 
620 
WRITE(*,*) (RMatrix(I,j),j=1,3) 
621 
END DO 
622 
END IF 
623  
624  
625 
! S matrix 
626 
S(1, 1) = Rmatrix(1, 1) + Rmatrix(2, 2) + Rmatrix(3, 3) 
627 
S(2, 1) = Rmatrix(2, 3)  Rmatrix(3, 2) 
628 
S(3, 1) = Rmatrix(3, 1)  Rmatrix(1, 3) 
629 
S(4, 1) = Rmatrix(1, 2)  Rmatrix(2, 1) 
630  
631 
S(1, 2) = S(2, 1) 
632 
S(2, 2) = Rmatrix(1, 1)  Rmatrix(2, 2)  Rmatrix(3, 3) 
633 
S(3, 2) = Rmatrix(1, 2) + Rmatrix(2, 1) 
634 
S(4, 2) = Rmatrix(1, 3) + Rmatrix(3, 1) 
635  
636 
S(1, 3) = S(3, 1) 
637 
S(2, 3) = S(3, 2) 
638 
S(3, 3) =Rmatrix(1, 1) + Rmatrix(2, 2)  Rmatrix(3, 3) 
639 
S(4, 3) = Rmatrix(2, 3) + Rmatrix(3, 2) 
640  
641 
S(1, 4) = S(4, 1) 
642 
S(2, 4) = S(4, 2) 
643 
S(3, 4) = S(4, 3) 
644 
S(4, 4) =Rmatrix(1, 1)  Rmatrix(2, 2) + Rmatrix(3, 3) 
645  
646  
647 
! PFL : I use my usual Jacobi diagonalisation 
648 
! Calculate eigenvalues and eigenvectors, and 
649 
! take the maximum eigenvalue lambda and the corresponding eigenvector q. 
650  
651 
IF (debug) THEN 
652 
WRITE(*,*) "S matrix" 
653 
DO I=1,4 
654 
WRITE(*,*) (S(I,j),j=1,4) 
655 
END DO 
656 
END IF 
657  
658 
Call Jacobi(S,4,EigVal,EigVec,4) 
659  
660 
Call Trie(4,EigVal,EigVec,4) 
661  
662 
Lambda=EigVal(4) 
663  
664 
! RMS Deviation 
665 
rmsd=sqrt(max(0.0d0,((x_norm+y_norm)2.0d0*lambda))/dble(na)) 
666  
667 
if (FRot.OR.FAlign) Call rotation_matrix(EigVec(1,4),U) 
668 
IF (FAlign) THEN 
669 
DO I=1,na 
670 
geom2(1,i)=Coord2(1,i)*U(1,1)+Coord2(2,i)*U(2,1)+Coord2(3,i)*U(3,1) +xc1 
671 
geom2(2,i)=Coord2(1,i)*U(1,2)+Coord2(2,i)*U(2,2)+Coord2(3,i)*U(3,2) +yc1 
672 
geom2(3,i)=Coord2(1,i)*U(1,3)+Coord2(2,i)*U(2,3)+Coord2(3,i)*U(3,3) +zc1 
673 
END DO 
674 
END IF 
675  
676 
END 
677  
678  
679 
! 
680 
subroutine rotation_matrix(q, U) 
681 
! 
682 
! This subroutine constructs rotation matrix U from quaternion q. 
683 
! 
684 
! This subroutine calculates RMSD using quaternions. 
685 
! It is based on the F90 routine bu E. Coutsias 
686 
! http://www.math.unm.edu/~vageli/homepage.html 
687 
! I (PFL) have just translated it, and I have changed the diagonalization 
688 
! subroutine. 
689 
! I also made some changes to make it suitable for Cart package. 
690 
! 
691 
! 
692 
! Copyright (C) 2004, 2005 Chaok Seok, Evangelos Coutsias and Ken Dill 
693 
! UCSF, Univeristy of New Mexico, Seoul National University 
694 
! Witten by Chaok Seok and Evangelos Coutsias 2004. 
695 

696 
! This library is free software; you can redistribute it and/or 
697 
! modify it under the terms of the GNU Lesser General Public 
698 
! License as published by the Free Software Foundation; either 
699 
! version 2.1 of the License, or (at your option) any later version. 
700 
! 
701 

702 
! This library is distributed in the hope that it will be useful, 
703 
! but WITHOUT ANY WARRANTY; without even the implied warranty of 
704 
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
705 
! Lesser General Public License for more details. 
706 
! 
707 

708 
! You should have received a copy of the GNU Lesser General Public 
709 
! License along with this library; if not, write to the Free Software 
710 
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA 
711 
! 
712  
713 
REAL(8) :: q(4) 
714 
REAL(8) :: U(3,3) 
715 
REAL(8) :: q0,q1,q2,q3,b0,b1,b2,b3,q00,q01,q02,q03 
716 
REAL(8) :: q11,q12,q13,q22,q23,q33 
717  
718 
q0 = q(1) 
719 
q1 = q(2) 
720 
q2 = q(3) 
721 
q3 = q(4) 
722  
723 
b0 = 2.0d0*q0 
724 
b1 = 2.0d0*q1 
725 
b2 = 2.0d0*q2 
726 
b3 = 2.0d0*q3 
727  
728 
q00 = b0*q01.0d0 
729 
q01 = b0*q1 
730 
q02 = b0*q2 
731 
q03 = b0*q3 
732  
733 
q11 = b1*q1 
734 
q12 = b1*q2 
735 
q13 = b1*q3 
736 

737 
q22 = b2*q2 
738 
q23 = b2*q3 
739 

740 
q33 = b3*q3 
741 

742 
U(1,1) = q00+q11 
743 
U(1,2) = q12q03 
744 
U(1,3) = q13+q02 
745 

746 
U(2,1) = q12+q03 
747 
U(2,2) = q00+q22 
748 
U(2,3) = q23q01 
749 

750 
U(3,1) = q13q02 
751 
U(3,2) = q23+q01 
752 
U(3,3) = q00+q33 
753 

754 
end 
755 

756 
!c============================================================ 
757 
!c 
758 
!c ++ diagonalisation par jacobi 
759 
!c vecteur propre i : V(i,i) 
760 
!c valeur propre i : D(i) 
761 
!c 
762 
!c============================================================ 
763 
!c 
764 
SUBROUTINE JACOBI(A,N,D,V,max_N) 
765 
IMPLICIT REAL(8) (AH,OZ) 
766 
parameter (max_it=500) 
767 
DIMENSION A(max_N,max_N),B(max_N),Z(max_N) 
768 
DIMENSION V(max_N,max_N),D(max_N) 
769  
770  
771 
DO 12 IP=1,N 
772 
DO 11 IQ=1,N 
773 
V(IP,IQ)=0. 
774 
11 CONTINUE 
775 
V(IP,IP)=1. 
776 
12 CONTINUE 
777 
DO 13 IP=1,N 
778 
B(IP)=A(IP,IP) 
779 
D(IP)=B(IP) 
780 
Z(IP)=0. 
781 
13 CONTINUE 
782 
NROT=0 
783 
DO 24 I=1,max_it 
784 
SM=0. 
785 
DO 15 IP=1,N1 
786 
DO 14 IQ=IP+1,N 
787 
SM=SM+ABS(A(IP,IQ)) 
788 
14 CONTINUE 
789 
15 CONTINUE 
790 
IF(SM.EQ.0.)RETURN 
791 
IF(I.LT.4)THEN 
792 
TRESH=0.2*SM/N**2 
793 
ELSE 
794 
TRESH=0. 
795 
ENDIF 
796 
DO 22 IP=1,N1 
797 
DO 21 IQ=IP+1,N 
798 
G=100.*ABS(A(IP,IQ)) 
799 
IF((I.GT.4).AND.(ABS(D(IP))+G.EQ.ABS(D(IP))) & 
800 
.AND.(ABS(D(IQ))+G.EQ.ABS(D(IQ))))THEN 
801 
A(IP,IQ)=0. 
802 
ELSE IF(ABS(A(IP,IQ)).GT.TRESH)THEN 
803 
H=D(IQ)D(IP) 
804 
IF(ABS(H)+G.EQ.ABS(H))THEN 
805 
T=A(IP,IQ)/H 
806 
ELSE 
807 
THETA=0.5*H/A(IP,IQ) 
808 
T=1./(ABS(THETA)+SQRT(1.+THETA**2)) 
809 
IF(THETA.LT.0.)T=T 
810 
ENDIF 
811 
C=1./SQRT(1+T**2) 
812 
S=T*C 
813 
TAU=S/(1.+C) 
814 
H=T*A(IP,IQ) 
815 
Z(IP)=Z(IP)H 
816 
Z(IQ)=Z(IQ)+H 
817 
D(IP)=D(IP)H 
818 
D(IQ)=D(IQ)+H 
819 
A(IP,IQ)=0. 
820 
DO 16 J=1,IP1 
821 
G=A(J,IP) 
822 
H=A(J,IQ) 
823 
A(J,IP)=GS*(H+G*TAU) 
824 
A(J,IQ)=H+S*(GH*TAU) 
825 
16 CONTINUE 
826 
DO 17 J=IP+1,IQ1 
827 
G=A(IP,J) 
828 
H=A(J,IQ) 
829 
A(IP,J)=GS*(H+G*TAU) 
830 
A(J,IQ)=H+S*(GH*TAU) 
831 
17 CONTINUE 
832 
DO 18 J=IQ+1,N 
833 
G=A(IP,J) 
834 
H=A(IQ,J) 
835 
A(IP,J)=GS*(H+G*TAU) 
836 
A(IQ,J)=H+S*(GH*TAU) 
837 
18 CONTINUE 
838 
DO 19 J=1,N 
839 
G=V(J,IP) 
840 
H=V(J,IQ) 
841 
V(J,IP)=GS*(H+G*TAU) 
842 
V(J,IQ)=H+S*(GH*TAU) 
843 
19 CONTINUE 
844 
NROT=NROT+1 
845 
ENDIF 
846 
21 CONTINUE 
847 
22 CONTINUE 
848 
DO 23 IP=1,N 
849 
B(IP)=B(IP)+Z(IP) 
850 
D(IP)=B(IP) 
851 
Z(IP)=0. 
852 
23 CONTINUE 
853 
24 CONTINUE 
854 
write(6,*) max_it,' iterations should never happen' 
855 
STOP 
856 
RETURN 
857 
END 
858 
! 
859 
!============================================================ 
860 
!c 
861 
!c ++ trie des vecteur dans l'ordre croissant 
862 
!c 
863 
!c============================================================ 
864 
!c 
865 
SUBROUTINE trie(nb_niv,ene,psi,max_niv) 
866 
integer(4) :: i,j,k,nb_niv,max_niv 
867 
REAL(8) :: ene(max_niv),psi(max_niv,max_niv) 
868 
REAL(8) :: a 
869  
870 
DO i=1,nb_niv 
871 
DO j=i+1,nb_niv 
872 
IF (ene(i) .GT. ene(j)) THEN 
873 
!c permutation 
874 
a=ene(i) 
875 
ene(i)=ene(j) 
876 
ene(j)=a 
877 
DO k=1,nb_niv 
878 
a=psi(k,i) 
879 
psi(k,i)=psi(k,j) 
880 
psi(k,j)=a 
881 
END DO 
882 
END IF 
883 
END DO 
884 
END DO 
885  
886 
END SUBROUTINE trie 
887 