root / utils / Xyz2Path.f90 @ 6
History | View | Annotate | Download (28.5 kB)
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 ',ng-1,' geometries' |
234 |
close(11) |
235 |
end |
236 |
|
237 |
!-------------------------------------------------------------- |
238 |
subroutine rtraitem(ifil,nnn,E,tab,atoms) |
239 |
! implicit REAL(8) :: (a-h,o-z) |
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:i-1) |
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*v2z-v1z*v2y |
401 |
v4y=v1z*v2x-v1x*v2z |
402 |
v4z=v1x*v2y-v1y*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*v2z-v5z*v2y) |
417 |
! * -v4y*(v5x*v2z-v5z*v2x) |
418 |
! * +v4z*(v5x*v2y-v5y*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(IC-32) |
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(IC-32) |
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=I-1 |
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 02110-1301 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 02110-1301 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*q0-1.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) = q12-q03 |
727 |
U(1,3) = q13+q02 |
728 |
|
729 |
U(2,1) = q12+q03 |
730 |
U(2,2) = q00+q22 |
731 |
U(2,3) = q23-q01 |
732 |
|
733 |
U(3,1) = q13-q02 |
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) (A-H,O-Z) |
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,N-1 |
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,N-1 |
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,IP-1 |
804 |
G=A(J,IP) |
805 |
H=A(J,IQ) |
806 |
A(J,IP)=G-S*(H+G*TAU) |
807 |
A(J,IQ)=H+S*(G-H*TAU) |
808 |
16 CONTINUE |
809 |
DO 17 J=IP+1,IQ-1 |
810 |
G=A(IP,J) |
811 |
H=A(J,IQ) |
812 |
A(IP,J)=G-S*(H+G*TAU) |
813 |
A(J,IQ)=H+S*(G-H*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)=G-S*(H+G*TAU) |
819 |
A(IQ,J)=H+S*(G-H*TAU) |
820 |
18 CONTINUE |
821 |
DO 19 J=1,N |
822 |
G=V(J,IP) |
823 |
H=V(J,IQ) |
824 |
V(J,IP)=G-S*(H+G*TAU) |
825 |
V(J,IQ)=H+S*(G-H*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 |
|