source: rmsdfun.f

Last change on this file was 38d77eb, checked in by baerbaer <baerbaer@…>, 14 years ago

Redirected standard out to logString.

SMMP produced a lot of log messages. This became an issue when run in massively
parallel environments. I replaced all writes to standard out to a write to logString.
The next step is to pass this string to a function that writes the messages to a log
file according to their importance and the chosen log level.

git-svn-id: svn+ssh://svn.berlios.de/svnroot/repos/smmp/trunk@34 26dc1dd8-5c4e-0410-9ffe-d298b4865968

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