1 | ! **************************************************************
|
---|
2 | !
|
---|
3 | ! This file contains the subroutines: rmsdfun,rmsdopt,fitmol,
|
---|
4 | ! jacobi,rmsinit
|
---|
5 | !
|
---|
6 | ! Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
7 | ! Shura Hayryan, Chin-Ku
|
---|
8 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
9 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
10 | !
|
---|
11 | ! **************************************************************
|
---|
12 |
|
---|
13 | real*8 function rmsdfun(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl)
|
---|
14 | !
|
---|
15 | ! --------------------------------------------------------------
|
---|
16 | ! Wrapping function for calculating rmsd
|
---|
17 | !
|
---|
18 | ! LIMITATION: requires call of rmsinit BEFORE calling this function
|
---|
19 | !
|
---|
20 | ! CALLS: rmsdopt
|
---|
21 | !
|
---|
22 | ! ---------------------------------------------------------------
|
---|
23 | !
|
---|
24 | include 'INCL.H'
|
---|
25 | include 'INCP.H'
|
---|
26 | !
|
---|
27 | ! Input
|
---|
28 | dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp)
|
---|
29 | ! Local
|
---|
30 | dimension rm(3,3),av1(3),av2(3)
|
---|
31 | call rmsdopt(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl,rm,av1,av2,rssd)
|
---|
32 |
|
---|
33 | rmsdfun = rssd
|
---|
34 |
|
---|
35 | return
|
---|
36 |
|
---|
37 | end
|
---|
38 |
|
---|
39 | !*******************************************************************
|
---|
40 | subroutine rmsdopt(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl,
|
---|
41 | & rm,av1,av2,rmsd)
|
---|
42 |
|
---|
43 | ! ---------------------------------------------------------------
|
---|
44 | ! PURPOSE: root mean square deviation (rmsd) between current SMMP
|
---|
45 | ! structure and reference atom coordinates 'x,y,zrf()'
|
---|
46 | ! for range of SMMP residues [ir1,ir2] in molecule 'nml'
|
---|
47 | !
|
---|
48 | ! ixat(i) - points to the atom in ref. coords., which is
|
---|
49 | ! equivalent to atom i of SMMP structure
|
---|
50 | ! (=0 if no equivalent in ref. structure exists)
|
---|
51 | !
|
---|
52 | ! isl = 0 : select all heavy atoms
|
---|
53 | ! isl = 1 : backbone atoms n,ca,c
|
---|
54 | ! isl = 2 : only ca atoms
|
---|
55 | !
|
---|
56 | ! CALLS: fitmol [S.K.Kearsley, Acta Cryst. 1989, A45, 208-210]
|
---|
57 | !
|
---|
58 | ! NB uncomment last lines in 'fitmol' to return coordinates
|
---|
59 | ! in 'x2' after fitting the ref. str. onto SMMP structure
|
---|
60 | ! ----------------------------------------------------------------
|
---|
61 |
|
---|
62 | include 'INCL.H'
|
---|
63 | include 'INCP.H'
|
---|
64 |
|
---|
65 | !-------------------------------------------------------- input
|
---|
66 | dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp)
|
---|
67 | !-------------------------------------------------------- output
|
---|
68 | dimension rm(3,3),av1(3),av2(3)
|
---|
69 | !-------------------------------------------------------- local
|
---|
70 | dimension x1(3,mxat),x2(3,mxat)
|
---|
71 | character*4 atnm
|
---|
72 |
|
---|
73 |
|
---|
74 | if (nml.lt.1.or.nml.gt.ntlml) then
|
---|
75 | write(*,*) ' rmsdopt> Sorry, there is no molecule #',nml
|
---|
76 | stop
|
---|
77 | endif
|
---|
78 |
|
---|
79 |
|
---|
80 | nr=0
|
---|
81 | na=0
|
---|
82 | n=0
|
---|
83 |
|
---|
84 | do im=1,ntlml
|
---|
85 | do ir=irsml1(im),irsml2(im)
|
---|
86 | if (im.eq.nml) nr=nr+1
|
---|
87 |
|
---|
88 | do ia=iatrs1(ir),iatrs2(ir)
|
---|
89 |
|
---|
90 | na=na+1
|
---|
91 |
|
---|
92 | if (im.eq.nml.and.nr.ge.ir1.and.nr.le.ir2) then ! range of res. for 'nml'
|
---|
93 |
|
---|
94 | ix=ixat(na)
|
---|
95 | atnm=nmat(ia)
|
---|
96 |
|
---|
97 | if (ix.gt.0.and.atnm(1:1).ne.'h') then
|
---|
98 |
|
---|
99 | if ( isl.eq.0
|
---|
100 |
|
---|
101 | & .or.
|
---|
102 |
|
---|
103 | & (isl.eq.1.and.(index(atnm,'n ').gt.0 .or.
|
---|
104 | & index(atnm,'ca ').gt.0 .or.
|
---|
105 | & index(atnm,'c ').gt.0 ))
|
---|
106 | & .or.
|
---|
107 |
|
---|
108 | & (isl.eq.2.and.index(atnm,'ca ').gt.0)
|
---|
109 |
|
---|
110 | & ) then
|
---|
111 |
|
---|
112 | n=n+1
|
---|
113 | x1(1,n)=xat(ia)
|
---|
114 | x1(2,n)=yat(ia)
|
---|
115 | x1(3,n)=zat(ia)
|
---|
116 | x2(1,n)=xrf(ix)
|
---|
117 | x2(2,n)=yrf(ix)
|
---|
118 | x2(3,n)=zrf(ix)
|
---|
119 |
|
---|
120 | endif ! atom selection
|
---|
121 | endif ! ix>0 & not 'h'
|
---|
122 |
|
---|
123 | endif ! res. range in mol. 'nml'
|
---|
124 |
|
---|
125 | enddo ! atoms
|
---|
126 | enddo ! residues
|
---|
127 | enddo ! molecules
|
---|
128 |
|
---|
129 | if (n.lt.3) then
|
---|
130 | write(*,*) ' rmsdopt> <3 atoms selected !'
|
---|
131 | stop
|
---|
132 | endif
|
---|
133 |
|
---|
134 | call fitmol(n,x1,x2, rm,av1,av2,rmsd)
|
---|
135 |
|
---|
136 | return
|
---|
137 | end
|
---|
138 | ! *********************************************
|
---|
139 | subroutine fitmol(n,x1,x2, rm,a1,a2,rmsd)
|
---|
140 | ! real*8 function fitmol(n,x1,x2)
|
---|
141 |
|
---|
142 | ! .......................................................
|
---|
143 | ! PURPOSE: compute RMSD of n positions in x1(3,) & x2(3,)
|
---|
144 | ! [S.K.Kearsley Acta Cryst. 1989,A45,208-210]
|
---|
145 | !
|
---|
146 | ! CALLS: jacobi
|
---|
147 | ! .......................................................
|
---|
148 | !f2py intent(out) rmsd
|
---|
149 |
|
---|
150 | include 'INCL.H'
|
---|
151 | ! implicit real*8 (a-h,o-z)
|
---|
152 | ! implicit integer*4 (i-n)
|
---|
153 |
|
---|
154 | ! ------------------------------------------- input/output
|
---|
155 | dimension x1(3,mxat),x2(3,mxat)
|
---|
156 | ! -------------------------------------------------- local
|
---|
157 | dimension e(4),q(4,4),v(4,4),dm(3),dp(3),a1(3),a2(3),rm(3,3)
|
---|
158 |
|
---|
159 | dn=dble(n)
|
---|
160 | ! ------------------- average of coordinates
|
---|
161 | do i=1,3
|
---|
162 | a1(i) = 0.d0
|
---|
163 | a2(i) = 0.d0
|
---|
164 | do j=1,n
|
---|
165 | a1(i) = a1(i) + x1(i,j)
|
---|
166 | a2(i) = a2(i) + x2(i,j)
|
---|
167 | enddo
|
---|
168 | a1(i) = a1(i)/dn
|
---|
169 | a2(i) = a2(i)/dn
|
---|
170 | enddo
|
---|
171 | ! ------------------------- compile quaternion
|
---|
172 | do i=1,4
|
---|
173 | do j=1,4
|
---|
174 | q(i,j)=0.d0
|
---|
175 | enddo
|
---|
176 | enddo
|
---|
177 |
|
---|
178 | do i=1,n
|
---|
179 |
|
---|
180 | do j=1,3
|
---|
181 | dm(j) = x1(j,i)-a1(j)
|
---|
182 | dp(j) = x2(j,i)-a2(j)
|
---|
183 | enddo
|
---|
184 |
|
---|
185 | dxm = dp(1) - dm(1)
|
---|
186 | dym = dp(2) - dm(2)
|
---|
187 | dzm = dp(3) - dm(3)
|
---|
188 | dxp = dp(1) + dm(1)
|
---|
189 | dyp = dp(2) + dm(2)
|
---|
190 | dzp = dp(3) + dm(3)
|
---|
191 |
|
---|
192 | q(1,1) = q(1,1) + dxm * dxm + dym * dym + dzm * dzm
|
---|
193 | q(1,2) = q(1,2) + dyp * dzm - dym * dzp
|
---|
194 | q(1,3) = q(1,3) + dxm * dzp - dxp * dzm
|
---|
195 | q(1,4) = q(1,4) + dxp * dym - dxm * dyp
|
---|
196 | q(2,2) = q(2,2) + dyp * dyp + dzp * dzp + dxm * dxm
|
---|
197 | q(2,3) = q(2,3) + dxm * dym - dxp * dyp
|
---|
198 | q(2,4) = q(2,4) + dxm * dzm - dxp * dzp
|
---|
199 | q(3,3) = q(3,3) + dxp * dxp + dzp * dzp + dym * dym
|
---|
200 | q(3,4) = q(3,4) + dym * dzm - dyp * dzp
|
---|
201 | q(4,4) = q(4,4) + dxp * dxp + dyp * dyp + dzm * dzm
|
---|
202 |
|
---|
203 | enddo
|
---|
204 |
|
---|
205 | do i=1,3
|
---|
206 | do j=i+1,4
|
---|
207 | q(j,i)=q(i,j)
|
---|
208 | enddo
|
---|
209 | enddo
|
---|
210 | ! ------------------------------ eigenvalues & -vectors
|
---|
211 | ndim4=4
|
---|
212 | call jacobi(q,ndim4,e,v)
|
---|
213 | ! --------------------------- lowest eigenvalue
|
---|
214 | im=1
|
---|
215 | em=e(1)
|
---|
216 | do i=2,4
|
---|
217 | if (e(i).lt.em) then
|
---|
218 | em=e(i)
|
---|
219 | im=i
|
---|
220 | endif
|
---|
221 | enddo
|
---|
222 |
|
---|
223 | rmsd = sqrt(em/dn)
|
---|
224 |
|
---|
225 | ! ================= uncomment following lines to fit molecule 2 onto 1
|
---|
226 |
|
---|
227 | ! ---------------------------------------------------rotation matrix
|
---|
228 | rm(1,1) = v(1,im)**2+v(2,im)**2-v(3,im)**2-v(4,im)**2
|
---|
229 | rm(1,2) = 2.d0*( v(2,im)*v(3,im)-v(1,im)*v(4,im) )
|
---|
230 | rm(1,3) = 2.d0*( v(2,im)*v(4,im)+v(1,im)*v(3,im) )
|
---|
231 | rm(2,1) = 2.d0*( v(2,im)*v(3,im)+v(1,im)*v(4,im) )
|
---|
232 | rm(2,2) = v(1,im)**2+v(3,im)**2-v(2,im)**2-v(4,im)**2
|
---|
233 | rm(2,3) = 2.d0*( v(3,im)*v(4,im)-v(1,im)*v(2,im) )
|
---|
234 | rm(3,1) = 2.d0*( v(2,im)*v(4,im)-v(1,im)*v(3,im) )
|
---|
235 | rm(3,2) = 2.d0*( v(3,im)*v(4,im)+v(1,im)*v(2,im) )
|
---|
236 | rm(3,3) = v(1,im)**2+v(4,im)**2-v(2,im)**2-v(3,im)**2
|
---|
237 |
|
---|
238 | ! do i=1,n
|
---|
239 | ! do j=1,3
|
---|
240 | ! dm(j) = x2(j,i) - a2(j)
|
---|
241 | ! enddo
|
---|
242 | ! do j=1,3
|
---|
243 | ! dp(j) = a1(j)
|
---|
244 | ! do k=1,3
|
---|
245 | ! dp(j) = dp(j) + rm(j,k) * dm(k)
|
---|
246 | ! enddo
|
---|
247 | ! x2(j,i) = dp(j)
|
---|
248 | ! enddo
|
---|
249 | ! enddo
|
---|
250 |
|
---|
251 | ! fitmol=rmsd
|
---|
252 |
|
---|
253 | return
|
---|
254 | end
|
---|
255 | ! ******************************
|
---|
256 | subroutine jacobi(a,n,d,v)
|
---|
257 |
|
---|
258 | ! ......................................................
|
---|
259 | ! PURPOSE: for given symmetric matrix 'a(n,n)
|
---|
260 | ! compute eigenvalues 'd' & eigenvectors 'v(,)'
|
---|
261 | !
|
---|
262 | ! [W.H.Press,S.A.Teukolsky,W.T.Vetterling,
|
---|
263 | ! B.P.Flannery, Numerical Recipes in FORTRAN,
|
---|
264 | ! Cambridge Univ. Press, 2nd Ed. 1992, 456-462]
|
---|
265 | !
|
---|
266 | ! CALLS: none
|
---|
267 | !
|
---|
268 | ! ......................................................
|
---|
269 |
|
---|
270 | !f2py intent(out) d
|
---|
271 | !f2py intent(out) v
|
---|
272 | parameter (NMAX=500)
|
---|
273 |
|
---|
274 | integer n,nrot,i,ip,iq,j
|
---|
275 |
|
---|
276 |
|
---|
277 | real*8 a(n,n),d(n),v(n,n),
|
---|
278 | & c,g,h,s,sm,t,tau,theta,tresh,b(NMAX),z(NMAX),smeps
|
---|
279 |
|
---|
280 | smeps=1.0d-6
|
---|
281 |
|
---|
282 |
|
---|
283 | do ip=1,n
|
---|
284 | b(ip)=a(ip,ip)
|
---|
285 | d(ip)=b(ip)
|
---|
286 | z(ip)=0.d0
|
---|
287 | do iq=1,n
|
---|
288 | v(ip,iq)=0.d0
|
---|
289 | enddo
|
---|
290 | v(ip,ip)=1.d0
|
---|
291 | enddo
|
---|
292 |
|
---|
293 | nrot=0
|
---|
294 |
|
---|
295 | do i=1,500
|
---|
296 |
|
---|
297 | sm=0.d0
|
---|
298 |
|
---|
299 | do ip=1,n-1
|
---|
300 | do iq=ip+1,n
|
---|
301 | sm=sm+abs(a(ip,iq))
|
---|
302 | enddo
|
---|
303 | enddo
|
---|
304 |
|
---|
305 | if (sm.le.smeps) return ! normal end
|
---|
306 |
|
---|
307 | if (i.lt.4) then
|
---|
308 | tresh=0.2d0*sm/n**2
|
---|
309 | else
|
---|
310 | tresh=0.d0
|
---|
311 | endif
|
---|
312 |
|
---|
313 | do ip=1,n-1
|
---|
314 | do iq=ip+1,n
|
---|
315 |
|
---|
316 | g=100.d0*abs(a(ip,iq))
|
---|
317 |
|
---|
318 | if((i.gt.4).and.(abs(d(ip))+
|
---|
319 |
|
---|
320 | &g.eq.abs(d(ip))).and.(abs(d(iq))+g.eq.abs(d(iq))))then
|
---|
321 | a(ip,iq)=0.d0
|
---|
322 |
|
---|
323 | else if(abs(a(ip,iq)).gt.tresh)then
|
---|
324 |
|
---|
325 | h=d(iq)-d(ip)
|
---|
326 |
|
---|
327 | if (abs(h)+g.eq.abs(h)) then
|
---|
328 |
|
---|
329 | t=a(ip,iq)/h
|
---|
330 |
|
---|
331 | else
|
---|
332 |
|
---|
333 | theta=0.5d0*h/a(ip,iq)
|
---|
334 | t=1.d0/(abs(theta)+sqrt(1.d0+theta**2))
|
---|
335 | if (theta.lt.0.d0) t=-t
|
---|
336 |
|
---|
337 | endif
|
---|
338 |
|
---|
339 | c=1.d0/sqrt(1.d0+t**2)
|
---|
340 | s=t*c
|
---|
341 | tau=s/(1.d0+c)
|
---|
342 | h=t*a(ip,iq)
|
---|
343 | z(ip)=z(ip)-h
|
---|
344 | z(iq)=z(iq)+h
|
---|
345 | d(ip)=d(ip)-h
|
---|
346 | d(iq)=d(iq)+h
|
---|
347 | a(ip,iq)=0.d0
|
---|
348 |
|
---|
349 | do j=1,ip-1
|
---|
350 | g=a(j,ip)
|
---|
351 | h=a(j,iq)
|
---|
352 | a(j,ip)=g-s*(h+g*tau)
|
---|
353 | a(j,iq)=h+s*(g-h*tau)
|
---|
354 | enddo
|
---|
355 |
|
---|
356 | do j=ip+1,iq-1
|
---|
357 | g=a(ip,j)
|
---|
358 | h=a(j,iq)
|
---|
359 | a(ip,j)=g-s*(h+g*tau)
|
---|
360 | a(j,iq)=h+s*(g-h*tau)
|
---|
361 | enddo
|
---|
362 | do j=iq+1,n
|
---|
363 | g=a(ip,j)
|
---|
364 | h=a(iq,j)
|
---|
365 | a(ip,j)=g-s*(h+g*tau)
|
---|
366 | a(iq,j)=h+s*(g-h*tau)
|
---|
367 | enddo
|
---|
368 | do j=1,n
|
---|
369 | g=v(j,ip)
|
---|
370 | h=v(j,iq)
|
---|
371 | v(j,ip)=g-s*(h+g*tau)
|
---|
372 | v(j,iq)=h+s*(g-h*tau)
|
---|
373 | enddo
|
---|
374 | nrot=nrot+1
|
---|
375 |
|
---|
376 | endif
|
---|
377 |
|
---|
378 | enddo
|
---|
379 | enddo
|
---|
380 |
|
---|
381 | do ip=1,n
|
---|
382 | b(ip)=b(ip)+z(ip)
|
---|
383 | d(ip)=b(ip)
|
---|
384 | z(ip)=0.d0
|
---|
385 | enddo
|
---|
386 |
|
---|
387 | enddo
|
---|
388 |
|
---|
389 | write(*,*) ' jacobi> too many iterations'
|
---|
390 | stop
|
---|
391 |
|
---|
392 | return
|
---|
393 | end
|
---|
394 |
|
---|
395 | ! ***********************************************************
|
---|
396 |
|
---|
397 | subroutine rmsinit(nml,string)
|
---|
398 | !
|
---|
399 | !------------------------------------------------------------------------------
|
---|
400 | ! Reads in pdb-file 'string' into INCP.H and initalizes
|
---|
401 | ! the files that 'rmdsopt' needs to calculate the rmsd
|
---|
402 | ! of a configuration with the pdb-configuration
|
---|
403 | !
|
---|
404 | ! CALLS: pdbread,atixpdb
|
---|
405 | !
|
---|
406 | ! ----------------------------------------------------------------------------
|
---|
407 | !
|
---|
408 | include 'INCL.H'
|
---|
409 | include 'INCP.H'
|
---|
410 |
|
---|
411 | character string*(*)
|
---|
412 |
|
---|
413 | if(string.eq.'smmp') then
|
---|
414 | !
|
---|
415 | ! Compare with a smmp-structure
|
---|
416 | !
|
---|
417 | do i=iatrs1(irsml1(nml)),iatrs2(irsml2(nml))
|
---|
418 | if(nmat(i)(1:1).ne.'h') then
|
---|
419 | ixatp(i)=i
|
---|
420 | else
|
---|
421 | ixatp(i) = 0
|
---|
422 | end if
|
---|
423 | enddo
|
---|
424 | !
|
---|
425 | else
|
---|
426 | !
|
---|
427 | ! Reference structure is read in from pdb-file
|
---|
428 | !
|
---|
429 | call pdbread(string,ier)
|
---|
430 | if(ier.ne.0) stop
|
---|
431 | call atixpdb(nml)
|
---|
432 | !
|
---|
433 | end if
|
---|
434 | print *,'RMSD initialized with ',string
|
---|
435 | return
|
---|
436 |
|
---|
437 | end
|
---|
438 |
|
---|