root / utils / Xyz2Path.f90 @ 6
Historique  Voir  Annoter  Télécharger (28,49 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 
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  
96 
AU2PS=1./41341.37 
97 
NbPrint=100 
98 
Pi=dacos(1.d0) 
99 
IOOUT=13 
100 
Conv=1. 
101 
FRot=.TRUE. 
102 
FAlign=.TRUE. 
103 
Debug=.FALSE. 
104 
NbDist=0 
105 
NbAngle=0 
106 
NbDie=0 
107  
108 
iarg=command_argument_count() 
109 
if (iarg.lt.1) then 
110 
write(*,*) 'XYZ filename:' 
111 
read(*,'(a)') f1 
112 
else 
113 
call getarg(1,f1) 
114 
endif 
115  
116 
open(13,file='Scan.dat') 
117  
118 
INQUIRE(File='list',Exist=FExist) 
119 
if (.NOT.FExist) THEN 
120 
WRITE(*,*) "No file 'list', just printing Energy" 
121 
END IF 
122  
123 
open(11,file=f1) 
124 
call rtraitem(11,nnn,ener,geos1,atoms) 
125 
close(11) 
126  
127 
DO I=1,nnn 
128 
Iat=ConvertNumAt(atoms(I)) 
129 
Mass(I)=MassAt(Iat) 
130 
! write(*,*) I,Atoms(I),Mass(I),geos1(:,I) 
131 
END DO 
132  
133 
Nat=nnn 
134  
135 
if (FExist) THEN 
136 
open(14,file='list') 
137 
Type="d" 
138 
DO WHILE (Type.ne."E") 
139 
CALL READLINE(14,Type,Line) 
140 
! WRITE(*,*) Line,Type 
141 
if (Type.eq."b") THEN 
142 
NbDist=NbDist+1 
143 
READ(Line,*) At1,At2 
144 
At1B(NbDist)=At1 
145 
At2B(NbDist)=At2 
146 
WRITE(13,'("# b ",2(A3,I3))') Atoms(At1),At1,Atoms(At2),At2 
147 
WRITE(*,'("# b ",2(A3,I3))') Atoms(At1),At1,Atoms(At2),At2 
148 
END IF 
149 
if (Type.eq."a") THEN 
150 
NbAngle=NbAngle+1 
151 
READ(Line,*) At1,At2,At3 
152 
At1A(NbAngle)=At1 
153 
At2A(NbAngle)=At2 
154 
At3A(NbAngle)=At3 
155 
WRITE(13,'("# a ",3(A3,I3))') Atoms(At1),At1,Atoms(At2), & 
156 
At2, Atoms(At3),At3 
157 
WRITE(*,'("# a ",3(A3,I3))') Atoms(At1),At1,Atoms(At2), & 
158 
At2, Atoms(At3),At3 
159 
END IF 
160 
if (Type.eq."d") THEN 
161 
NbDie=NbDie+1 
162 
READ(Line,*) At1,At2,At3,At4 
163 
At1D(NbDie)=At1 
164 
At2D(NbDie)=At2 
165 
At3D(NbDie)=At3 
166 
At4D(NbDie)=At4 
167 
WRITE(13,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2, & 
168 
Atoms(At3),At3,Atoms(At4),At4 
169 
WRITE(*,'("# d ",4(A3,I3))') Atoms(At1),At1,Atoms(At2),At2, & 
170 
Atoms(At3),At3,Atoms(At4),At4 
171  
172 
END IF 
173 
if (Type.eq."c") THEN 
174 
NbCOM=NbCOM+1 
175 
READ(Line,*) At1,(AtCom(j,NbCOM),j=1,At1) 
176 
AtCom(0,NbCOM)=At1 
177 
Atoms(nat+NbCom)="G" 
178 
WRITE(13,'("# c ",I4,20(A3,I3))') At1, & 
179 
(Atoms(AtCom(i,NbCoM)),AtCom(i,NbCOM),i=1,At1) 
180 
WRITE(*,'("# c ",I4,20(A3,I3))') At1, & 
181 
(Atoms(AtCom(i,NbCoM)),AtCom(i,NbCOM),i=1,At1) 
182 
END IF 
183  
184 
END DO 
185 
END IF 
186  
187  
188  
189 
fmt='( (1X,F12.5),1X,F15.6)' 
190 
write(fmt(2:4),'(i3)') NbDist+NbAngle+NbDie+1 
191 
! write(*,*) nat3 
192 
! write(*,*) 'fmt:',fmt 
193  
194 
ng=1 
195  
196 
s=0. 
197  
198 
open(11,file=f1) 
199  
200 
10 call rtraitem(11,nnn,ener,geos,atoms) 
201 
! WRITE(*,*) nnn 
202 
if (nnn.gt.0) then 
203  
204 
call CalcRmsd(nnn, geos1, geos,MRot,rmsd,FRot,FAlign,debug) 
205 
ds=0. 
206 
DO I=1,nnn 
207 
DO J=1,3 
208 
ds=ds+mass(I)*(geos(J,I)geos1(J,I))**2 
209 
! write(*,*) I,J,geos(J,I),Geos1(J,I),ds 
210 
geos1(J,I)=Geos(J,I) 
211 
END DO 
212 
END DO 
213 
s=s+sqrt(ds) 
214  
215 
! We convert coordinates in a0 into Angstroem 
216 
! write(*,*) "Analyse !" 
217 
if (FExist) THEN 
218 
Call Analyse(geos) 
219  
220 
write(IOOUT,fmt) s,(VB(j)/Conv,j=1,NbDist), & 
221 
(VA(j),j=1,NbAngle), & 
222 
(VD(j),j=1,NbDie),ener 
223 
write(*,fmt) s,(VB(j)/Conv,j=1,NbDist), & 
224 
(VA(j),j=1,NbAngle), & 
225 
(VD(j),j=1,NbDie),ener 
226 
ELSE 
227 
write(IOOUT,fmt) s,ener 
228 
write(*,fmt) s,ener 
229 
END IF 
230 
ng=ng+1 
231 
goto 10 
232 
endif 
233 
WRITE(*,*) 'Found ',ng1,' geometries' 
234 
close(11) 
235 
end 
236  
237 
! 
238 
subroutine rtraitem(ifil,nnn,E,tab,atoms) 
239 
! implicit REAL(8) :: (ah,oz) 
240 
IMPLICIT NONE 
241  
242 
character(120) :: line 
243 
character(3) :: Atoms(*) 
244 
INTEGER(4) :: nnn, IFil, Idx, I, J 
245 
REAL(8) :: Tab(3,*), E 
246  
247  
248 
read(ifil,*,err=99,end=99) nnn 
249 
read(ifil,'(A)') line 
250 
idx=index(line,'E') 
251 
! WRITE(*,*) 'idx,line',idx,line 
252 
if (idx/=0) THEN 
253 
line=line(idx+2:) 
254 
idx=index(line,"=") 
255 
if (idx/=0) line=line(idx+1:) 
256 
idx=index(line,":") 
257 
if (idx/=0) line=line(idx+1:) 
258 
read(line,*) E 
259 
ELSE 
260 
E=0. 
261 
END IF 
262 

263 
! write(*,*) 'coucou',line 
264 
do i=1,nnn 
265 
read(ifil,'(A)',err=99,end=99) Line 
266 
do while (line(1:1)==' ') 
267 
line=line(2:) 
268 
end do 
269 
atoms(i)=line(1:3) 
270 
! write(*,*) 'coucou atoms',atoms(i) 
271 
do while (line(1:1).ne.' ') 
272 
line=line(2:) 
273 
end do 
274 
! WRITE(*,*) 'coucou2:',i,line 
275 
read(line,*) (tab(j,i),j=1,3) 
276 
end do 
277 
! WRITE(*,*) 'coucou:',nnn,tab(1,1) 
278 
return 
279 
99 nnn=0 
280 
! WRITE(*,*) 'Erreur lecture',ifil,nnn 
281 
return 
282 
end 
283  
284 
! 
285  
286 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
287 
! READLINE 
288 
! This subroutine reads a line for a file, and converts 
289 
! the first field into a character variable, and the rest into 4 integers 
290 
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
291  
292 
SUBROUTINE READLINE(IOIN,Type,Line) 
293 

294 
IMPLICIT NONE 
295 
CHARACTER(5) :: Type 
296 
CHARACTER(120) :: LINE 
297 
INTEGER(4) :: i,IOIN 
298 
READ(IOIN,'(A120)',ERR=999,END=999) LINE 
299 
! WRITE(*,*) Line 
300 
DO WHILE (LINE(1:1).eq.' ') 
301 
LINE=LINE(2:) 
302 
END DO 
303  
304 
i=1 
305 
DO WHILE (LINE(i:i).ne.' ') 
306 
i=i+1 
307 
END DO 
308 

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

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

494 
! This library is free software; you can redistribute it and/or 
495 
! modify it under the terms of the GNU Lesser General Public 
496 
! License as published by the Free Software Foundation; either 
497 
! version 2.1 of the License, or (at your option) any later version. 
498 
! 
499 

500 
! This library is distributed in the hope that it will be useful, 
501 
! but WITHOUT ANY WARRANTY; without even the implied warranty of 
502 
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
503 
! Lesser General Public License for more details. 
504 
! 
505 

506 
! You should have received a copy of the GNU Lesser General Public 
507 
! License along with this library; if not, write to the Free Software 
508 
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA 
509 
! 
510  
511 
subroutine CalcRmsd(na, geom,geom2,U,rmsd,FRot,FAlign,debug) 
512 
! 
513 
! This subroutine calculates the least square rmsd of two coordinate 
514 
! sets coord1(3,n) and coord2(3,n) using a method based on quaternion. 
515 
! It then calculate the rotation matrix U and the centers of coord, and uses 
516 
! them to align the two molecules. 
517 
! 
518  
519  
520 
IMPLICIT NONE 
521  
522 
INCLUDE "Xyz2Path.param" 
523  
524 
INTEGER(4) :: IOIN, IOOUT, IOSCRN, IOPRNT, IOTAMP 
525 
COMMON/IODEFS/IOIN,IOOUT,IOSCRN,IOPRNT,IOTAMP 
526  
527 
INTEGER(4) :: na 
528 
REAL(8) :: geom(3,MaxNAt), Geom2(3,MaxNAt) 
529 
REAL(8) :: U(3,3), rmsd 
530 
LOGICAL FRot,FAlign,Debug 
531  
532 
REAL(8) :: Coord1(3,MaxNAt), Coord2(3,MaxNAt) 
533 
REAL(8) :: xc1,yc1,zc1, xc2,yc2,zc2 
534 

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

679 
! This library is free software; you can redistribute it and/or 
680 
! modify it under the terms of the GNU Lesser General Public 
681 
! License as published by the Free Software Foundation; either 
682 
! version 2.1 of the License, or (at your option) any later version. 
683 
! 
684 

685 
! This library is distributed in the hope that it will be useful, 
686 
! but WITHOUT ANY WARRANTY; without even the implied warranty of 
687 
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 
688 
! Lesser General Public License for more details. 
689 
! 
690 

691 
! You should have received a copy of the GNU Lesser General Public 
692 
! License along with this library; if not, write to the Free Software 
693 
! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 021101301 USA 
694 
! 
695  
696 
REAL(8) :: q(4) 
697 
REAL(8) :: U(3,3) 
698 
REAL(8) :: q0,q1,q2,q3,b0,b1,b2,b3,q00,q01,q02,q03 
699 
REAL(8) :: q11,q12,q13,q22,q23,q33 
700  
701 
q0 = q(1) 
702 
q1 = q(2) 
703 
q2 = q(3) 
704 
q3 = q(4) 
705  
706 
b0 = 2.0d0*q0 
707 
b1 = 2.0d0*q1 
708 
b2 = 2.0d0*q2 
709 
b3 = 2.0d0*q3 
710  
711 
q00 = b0*q01.0d0 
712 
q01 = b0*q1 
713 
q02 = b0*q2 
714 
q03 = b0*q3 
715  
716 
q11 = b1*q1 
717 
q12 = b1*q2 
718 
q13 = b1*q3 
719 

720 
q22 = b2*q2 
721 
q23 = b2*q3 
722 

723 
q33 = b3*q3 
724 

725 
U(1,1) = q00+q11 
726 
U(1,2) = q12q03 
727 
U(1,3) = q13+q02 
728 

729 
U(2,1) = q12+q03 
730 
U(2,2) = q00+q22 
731 
U(2,3) = q23q01 
732 

733 
U(3,1) = q13q02 
734 
U(3,2) = q23+q01 
735 
U(3,3) = q00+q33 
736 

737 
end 
738 

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