[e40e335] | 1 | c **************************************************************
|
---|
| 2 | c
|
---|
| 3 | c This file contains the subroutines: rmsdfun,rmsdopt,fitmol,
|
---|
| 4 | c jacobi,rmsinit
|
---|
| 5 | c
|
---|
| 6 | c Copyright 2003-2005 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
| 7 | c Shura Hayryan, Chin-Ku
|
---|
| 8 | c Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
| 9 | c Jan H. Meinke, Sandipan Mohanty
|
---|
| 10 | c
|
---|
| 11 | c **************************************************************
|
---|
| 12 |
|
---|
| 13 | real*8 function rmsdfun(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl)
|
---|
| 14 | C
|
---|
| 15 | C --------------------------------------------------------------
|
---|
| 16 | c Wrapping function for calculating rmsd
|
---|
| 17 | c
|
---|
| 18 | c LIMITATION: requires call of rmsinit BEFORE calling this function
|
---|
| 19 | c
|
---|
| 20 | c CALLS: rmsdopt
|
---|
| 21 | c
|
---|
| 22 | c ---------------------------------------------------------------
|
---|
| 23 | c
|
---|
| 24 | include 'INCL.H'
|
---|
| 25 | include 'INCP.H'
|
---|
| 26 | c
|
---|
| 27 | c Input
|
---|
| 28 | dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp)
|
---|
| 29 | c 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 | c*******************************************************************
|
---|
| 40 | subroutine rmsdopt(nml,ir1,ir2,ixat,xrf,yrf,zrf,isl,
|
---|
| 41 | # rm,av1,av2,rmsd)
|
---|
| 42 |
|
---|
| 43 | c ---------------------------------------------------------------
|
---|
| 44 | c PURPOSE: root mean square deviation (rmsd) between current SMMP
|
---|
| 45 | c structure and reference atom coordinates 'x,y,zrf()'
|
---|
| 46 | c for range of SMMP residues [ir1,ir2] in molecule 'nml'
|
---|
| 47 | c
|
---|
| 48 | c ixat(i) - points to the atom in ref. coords., which is
|
---|
| 49 | c equivalent to atom i of SMMP structure
|
---|
| 50 | c (=0 if no equivalent in ref. structure exists)
|
---|
| 51 | c
|
---|
| 52 | c isl = 0 : select all heavy atoms
|
---|
| 53 | c isl = 1 : backbone atoms n,ca,c
|
---|
| 54 | c isl = 2 : only ca atoms
|
---|
| 55 | c
|
---|
| 56 | c CALLS: fitmol [S.K.Kearsley, Acta Cryst. 1989, A45, 208-210]
|
---|
| 57 | c
|
---|
| 58 | c NB uncomment last lines in 'fitmol' to return coordinates
|
---|
| 59 | c in 'x2' after fitting the ref. str. onto SMMP structure
|
---|
| 60 | c ----------------------------------------------------------------
|
---|
| 61 |
|
---|
| 62 | include 'INCL.H'
|
---|
| 63 | include 'INCP.H'
|
---|
| 64 |
|
---|
| 65 | c-------------------------------------------------------- input
|
---|
| 66 | dimension ixat(mxat),xrf(mxatp),yrf(mxatp),zrf(mxatp)
|
---|
| 67 | c-------------------------------------------------------- output
|
---|
| 68 | dimension rm(3,3),av1(3),av2(3)
|
---|
| 69 | c-------------------------------------------------------- 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 | c *********************************************
|
---|
| 139 | subroutine fitmol(n,x1,x2, rm,a1,a2,rmsd)
|
---|
| 140 | c real*8 function fitmol(n,x1,x2)
|
---|
| 141 |
|
---|
| 142 | c .......................................................
|
---|
| 143 | c PURPOSE: compute RMSD of n positions in x1(3,) & x2(3,)
|
---|
| 144 | c [S.K.Kearsley Acta Cryst. 1989,A45,208-210]
|
---|
| 145 | c
|
---|
| 146 | c CALLS: jacobi
|
---|
| 147 | c .......................................................
|
---|
| 148 | cf2py intent(out) rmsd
|
---|
| 149 |
|
---|
| 150 | include 'INCL.H'
|
---|
| 151 | c implicit real*8 (a-h,o-z)
|
---|
| 152 | c implicit integer*4 (i-n)
|
---|
| 153 |
|
---|
| 154 | c ------------------------------------------- input/output
|
---|
| 155 | dimension x1(3,mxat),x2(3,mxat)
|
---|
| 156 | c -------------------------------------------------- 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 | c ------------------- 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 | c ------------------------- 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 | c ------------------------------ eigenvalues & -vectors
|
---|
| 211 | ndim4=4
|
---|
| 212 | call jacobi(q,ndim4,e,v)
|
---|
| 213 | c --------------------------- 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 | c ================= uncomment following lines to fit molecule 2 onto 1
|
---|
| 226 |
|
---|
| 227 | c ---------------------------------------------------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 | c do i=1,n
|
---|
| 239 | c do j=1,3
|
---|
| 240 | c dm(j) = x2(j,i) - a2(j)
|
---|
| 241 | c enddo
|
---|
| 242 | c do j=1,3
|
---|
| 243 | c dp(j) = a1(j)
|
---|
| 244 | c do k=1,3
|
---|
| 245 | c dp(j) = dp(j) + rm(j,k) * dm(k)
|
---|
| 246 | c enddo
|
---|
| 247 | c x2(j,i) = dp(j)
|
---|
| 248 | c enddo
|
---|
| 249 | c enddo
|
---|
| 250 |
|
---|
| 251 | c fitmol=rmsd
|
---|
| 252 |
|
---|
| 253 | return
|
---|
| 254 | end
|
---|
| 255 | c ******************************
|
---|
| 256 | subroutine jacobi(a,n,d,v)
|
---|
| 257 |
|
---|
| 258 | c ......................................................
|
---|
| 259 | c PURPOSE: for given symmetric matrix 'a(n,n)
|
---|
| 260 | c compute eigenvalues 'd' & eigenvectors 'v(,)'
|
---|
| 261 | c
|
---|
| 262 | c [W.H.Press,S.A.Teukolsky,W.T.Vetterling,
|
---|
| 263 | c B.P.Flannery, Numerical Recipes in FORTRAN,
|
---|
| 264 | c Cambridge Univ. Press, 2nd Ed. 1992, 456-462]
|
---|
| 265 | c
|
---|
| 266 | c CALLS: none
|
---|
| 267 | c
|
---|
| 268 | c ......................................................
|
---|
| 269 |
|
---|
| 270 | cf2py intent(out) d
|
---|
| 271 | cf2py 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 | c ***********************************************************
|
---|
| 396 |
|
---|
| 397 | subroutine rmsinit(nml,string)
|
---|
| 398 | c
|
---|
| 399 | c------------------------------------------------------------------------------
|
---|
| 400 | c Reads in pdb-file 'string' into INCP.H and initalizes
|
---|
| 401 | c the files that 'rmdsopt' needs to calculate the rmsd
|
---|
| 402 | c of a configuration with the pdb-configuration
|
---|
| 403 | C
|
---|
| 404 | c CALLS: pdbread,atixpdb
|
---|
| 405 | c
|
---|
| 406 | c ----------------------------------------------------------------------------
|
---|
| 407 | c
|
---|
| 408 | include 'INCL.H'
|
---|
| 409 | include 'INCP.H'
|
---|
| 410 |
|
---|
| 411 | character string*(*)
|
---|
| 412 |
|
---|
| 413 | if(string.eq.'smmp') then
|
---|
| 414 | c
|
---|
| 415 | c Compare with a smmp-structure
|
---|
| 416 | c
|
---|
| 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 | c
|
---|
| 425 | else
|
---|
| 426 | c
|
---|
| 427 | c Reference structure is read in from pdb-file
|
---|
| 428 | c
|
---|
| 429 | call pdbread(string,ier)
|
---|
| 430 | if(ier.ne.0) stop
|
---|
| 431 | call atixpdb(nml)
|
---|
| 432 | c
|
---|
| 433 | end if
|
---|
| 434 | print *,'RMSD initialized with ',string
|
---|
| 435 | return
|
---|
| 436 |
|
---|
| 437 | end
|
---|
| 438 |
|
---|