| 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 |
|
|---|
