source/fragments/suppdb_ca_dis.f

   1       dimension s1(3,10000),s2(3,10000),ca1(3,10000),ca2(3,10000)
   2       dimension xyz1(3,5000),xyz2(3,5000),waby(5000),ires(5000),
   3      & atnam(5000),rms(500,500),imin(5000),imax(5000),isel(5000),
   4      & isel_t(5000),issel(5000),diff(5000),iperm(5000)
   5       character*4 atnam
   6       dimension t(3),b(3,3),z(3,5000)
   7       integer resnum,resnum1,resnum2,resnum3,resnum4,frag(5000)
   8       character*80 karta
   9       character*80 fnam1,fnam2,fnam3
  10       logical insel
  11       ifrag=1
  12       nnn = iargc()
  13       if (nnn.lt.3) then
  14         stop
  15      & 'Usage: suppdb pdb_file_1 pdb_file_2 out_file [resnum1 resnum2]'
  16       endif
  17       call getarg(1,fnam1)
  18       call getarg(2,fnam2)
  19       call getarg(3,fnam3)
  20       open (2,file=fnam3,status='unknown')
  21       if (nnn.lt.5) then
  22         print *,'resnum1,resnum2'
  23         read (*,*) resnum1,resnum2
  24       else
  25         call getarg(4,karta)
  26         read (karta,*) resnum1
  27         call getarg(5,karta)
  28         read (karta,*) resnum2
  29       endif
  30       if (nnn.gt.5) then
  31         call getarg(6,karta)
  32         read (karta,*) resnum3
  33       else
  34         resnum3=resnum1
  35       endif
  36       resnum4=resnum3+resnum2-resnum1
  37       print *,'resnum1',resnum1,' resnum2',resnum2
  38       print *,'resnum3',resnum3,' resnum4',resnum4
  39       print *,'opening file ',fnam1
  40       open (1,file=fnam1,status='old',err=111)
  41       iat=0
  42       isup=0
  43       write (2,'(/2a)') 'Coordinates of the first set from ',
  44      &  fnam1(:len_trim(fnam1))
  45       do irec=1,10000
  46       read (1,'(a80)',end=10) karta
  47       if ((karta(1:4).eq.'ATOM'.or.karta(1:6).eq.'HETATM').and.
  48      & karta(14:15).eq.'CA') then
  49       read (karta(23:26),*) resnum
  50       if (resnum.ge.resnum1 .and. resnum.le.resnum2) then
  51       iat=iat+1
  52       read (karta(14:16),'(a3)') atnam(iat)
  53       read (karta(31:54),'(3f8.3)') xyz1(1,iat),xyz1(2,iat),
  54      1                              xyz1(3,iat)
  55       read (karta(23:26),'(i4)') ires(iat)
  56       read (karta(18:20),'(a3)') waby(iat)
  57       if (atnam(iat).eq.'CA ') then
  58         write (2,'(a)') karta
  59         isup=isup+1
  60         do i=1,3
  61           ca1(i,isup)=xyz1(i,iat)
  62           s1(i,isup)=xyz1(i,iat)
  63         enddo
  64       endif
  65       endif
  66       endif
  67       enddo
  68    10 nsup=isup
  69       nat=iat
  70       print *,'from first file: nat=',nat
  71       close(1)
  72       print *,'opening file ',fnam2
  73       open (1,file=fnam2,status='old',err=112)
  74       iat=0
  75       isup=0
  76       write (2,'(/2a)') 'Coordinates of the first set from ',
  77      &  fnam2(:len_trim(fnam2))
  78       do irec=1,10000
  79       read (1,'(a80)',end=20) karta
  80       if ((karta(1:4).eq.'ATOM'.or.karta(1:6).eq.'HETATM').and.
  81      & karta(14:15).eq.'CA') then
  82       read (karta(23:26),*) resnum
  83       if (resnum.ge.resnum3 .and. resnum.le.resnum4) then
  84       iat=iat+1
  85       read (karta(31:54),'(3f8.3)') xyz2(1,iat),xyz2(2,iat),
  86      1                              xyz2(3,iat)
  87       if (atnam(iat).eq.'CA ') then
  88         write (2,'(a)') karta
  89         isup=isup+1
  90         do i=1,3
  91           ca2(i,isup)=xyz2(i,iat)
  92           s2(i,isup)=xyz2(i,iat)
  93         enddo
  94       endif
  95       endif
  96       endif
  97       enddo
  98    20 continue
  99       print *,'from second file: nat=',iat
 100 !      write (2,'(//4a8)') '# length','  start','   rms'
 101
 102
 103       do i=1,nsup
 104         issel(i)=i
 105       enddo
 106
 107       nsup_t = nsup
 108
 109       nsup_old=0
 110
 111       DO WHILE (NSUP>4 .and. nsup.ne.nsup_old)
 112
 113       nsup_old = nsup
 114
 115       do is1=4,nsup
 116 !        write (2,'(a)')
 117         do is2=1,nsup-is1+1
 118           call fitsq(rms(is1,is2),s2(1,is2),s1(1,is2),is1,t,b)
 119 !          write (2,'(2(i4,4x),2f8.2)') is1,is2+resnum1-1,rms(is1,is2)
 120         enddo
 121       enddo
 122       write (2,'(//2a8,3a8)') '# length',' start','rmsmin',
 123      &  'rmsmax','rmsave'
 124       do i=4,nsup
 125         imin(i)=1
 126         imax(i)=1
 127         rmin=rms(i,1)
 128         rmax=rms(i,1)
 129         rmsave=rms(i,1)**2
 130         do j=2,nsup-i+1
 131           if ( rms(i,j).lt.rmin) then
 132             imin(i)=j
 133             rmin=rms(i,j)
 134           else if ( rms(i,j).gt.rmax) then
 135             imax(i)=j
 136             rmax=rms(i,j)
 137           endif
 138           rmsave=rmsave+rms(i,j)*rms(i,j)
 139         enddo
 140         rmsave=sqrt(rmsave/(nsup-i+1))
 141         write (2,'(2(i4,4x),3(f8.2))') i,imin(i)+resnum1-1,
 142      &   rms(i,imin(i)), ! rms(i,imin(i))/sqrt(float(i)),
 143      &   rms(i,imax(i)), ! rms(i,imax(i))/sqrt(float(i)),
 144      &   rmsave ! ,rmsave/sqrt(float(i))
 145       enddo
 146       rmscut=7.0d0
 147       rmscut1=1.0d0*rmscut
 148       do i4=nsup,4,-1
 149         if (rms(i4,imin(i4)).le.rmscut) exit
 150       enddo
 151       call fitsq(rrms,s2(1,imin(i4)),s1(1,imin(i4)),i4,t,b)
 152       nsup_t=nsup
 153       do i=1,nsup
 154         isel(i)=issel(i)
 155       enddo
 156       rrms=1000.0d0
 157 c      do while (rrms.gt.rmscut)
 158       do while (rrms.gt.rmscut)
 159         nnsup_t=0
 160         std=0.0d0
 161         diffmax=-10000.0d0
 162         do i=1,nsup_t
 163           ii=isel(i)
 164           call mvvad(b,ca2(1,ii),t,z(1,i))
 165           diffi=0.0d0
 166           do j=1,3
 167             diffi=diffi+(z(j,i)-ca1(j,ii))**2
 168           enddo
 169           std=std+diffi
 170           diffi=sqrt(diffi)
 171 !          write (2,*) i,ii,diffi
 172           diff(i)=diffi
 173           if (diffi.gt.diffmax) diffmax=diffi
 174         enddo
 175         do i=1,nsup_t
 176           iperm(i)=i
 177         enddo
 178         write(2,*) 'diff'
 179         write(2,*) (diff(i),i=1,nsup_t)
 180         call mysort(nsup_t,diff,iperm)
 181         write(2,*) 'diff'
 182         write(2,*) (diff(i),i=1,nsup_t)
 183         std=sqrt(std/nsup_t)
 184         write (2,*) 'std=',std
 185         write (*,*) 'std=',std
 186         sstd=dmin1(1.5d0*std,0.9999d0*diffmax)
 187 !        do i=1,nsup_t
 188 !          ii=isel(i)
 189 !          if (diff(i).lt.sstd .or. diff(i).lt.rmscut) then
 190 !            nnsup_t=nnsup_t+1
 191 !            isel(nnsup_t)=ii
 192 !          endif
 193 !        enddo
 194         sumdiff=0.0d0
 195         do i=1,nsup_t
 196           ii=isel(iperm(i))
 197 c          sumdiff=sumdiff+diff(i)**2
 198 c          if (dsqrt(sumdiff/dfloat(i)).lt.rmscut1) then
 199           if (diff(i).lt.rmscut1) then
 200             nnsup_t=nnsup_t+1
 201             isel_t(nnsup_t)=ii
 202           endif
 203         enddo
 204         nsup_t=nnsup_t
 205         do i=1,nsup_t
 206           iperm(i)=i
 207           isel(i)=isel_t(i)
 208         enddo
 209         if (nsup_t.eq.0 .or. nnsup_t.eq.0) exit
 210         write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
 211         write (2,'(20i4)') (isel(i),i=1,nnsup_t)
 212         call mysort(nsup_t,isel,iperm)
 213         print *,'nsup_t',nsup_t,' nnsup_t',nnsup_t
 214         print '(20i4)',(isel(i),i=1,nnsup_t)
 215         write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
 216         write (2,'(20i4)') (isel(i),i=1,nnsup_t)
 217         len_min=4
 218         i=1
 219         len=1
 220         do i=1,nsup_t
 221           ii=isel(i)
 222           do j=1,3
 223             s1(j,i)=ca1(j,ii)
 224             s2(j,i)=ca2(j,ii)
 225           enddo
 226         enddo
 227         do i=1,nsup_t
 228           ii=isel(i)
 229           do j=1,3
 230             s1(j,i)=ca1(j,ii)
 231             s2(j,i)=ca2(j,ii)
 232           enddo
 233         enddo
 234         if (nsup_t.le.2) exit
 235         print *,"before fitsq: nsup_r",nsup_t
 236         call fitsq(rrms,s2(1,1),s1(1,1),nsup_t,t,b)
 237         write (2,*) 'nsup_t=',nsup_t,' rrms=',rrms
 238         write (*,*) 'nsup_t=',nsup_t,' rrms=',rrms
 239       enddo
 240       i=1
 241       do while (i.le.nnsup_t)
 242         write (*,*) 'i=',i,' isel',isel(i),' isel1',isel(i+1)
 243         if (i.eq.nnsup_t .or. isel(i+1)-isel(i).gt.1) then
 244           if (len.lt.len_min) then
 245             write (*,'(20i4)') (isel(k),k=1,nnsup_t)
 246             write (*,*) 'len=',len,' len_min=',len_min
 247             do j=i+1,nnsup_t
 248               isel(j-len)=isel(j)
 249             enddo
 250             nnsup_t=nnsup_t-len
 251 !            write (2,'(20i4)') (isel(k),k=1,nnsup_t)
 252             i=i-len+1
 253             len=1
 254           else
 255             len=1
 256             i=i+1
 257           endif
 258         else
 259           len=len+1
 260           i=i+1
 261         endif
 262 !        write (2,*) 'i=',i,' len=',len
 263       enddo
 264       print *,'final nsup_t',nsup_t
 265       print *,"final isel"
 266       print '(20i4)',(isel(i),i=1,nnsup_t)
 267       write (2,*) 'final nsup_t',nsup_t
 268       write (2,*) "final isel"
 269       write (2,'(20i4)') (isel(i),i=1,nnsup_t)
 270       if (nnsup_t.gt.10) then
 271        do i=1,nnsup_t
 272         frag(isel(i))=ifrag
 273        enddo
 274        ifrag=ifrag+1
 275       endif
 276       print *,"----- end of iteration"
 277       write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
 278       nsup_t=nnsup_t
 279       write (2,'(/4a)') 'REMARK superimposed coordinates from files ',
 280      &fnam1(:len_trim(fnam1)),' and ',fnam2(:len_trim(fnam2))
 281       write (2,'(a,f6.2,a)')
 282      & 'REMARK Calpha rms deviation:',rrms,
 283      & ' angstroms'
 284       write (2,'(a,I5)') 'REMARK number of superposed residues',nsup_t
 285       i1=isel(1)
 286       do i=2,nnsup_t
 287         if (isel(i)-isel(i-1).gt.1) then
 288           write (2,'(a,i4,a,i4,a,f5.1,a)')
 289      &     'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1
 290            i1=isel(i)
 291         endif
 292       enddo
 293       write (2,'(a,i4,a,i4,a,f5.1,a)')
 294      &     'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1
 295       do i=1,nsup_t
 296         ii=isel(i)
 297         write (2,100) ii+resnum1-1,atnam(ii),waby(ii),ires(ii),
 298      &                (xyz1(j,ii),j=1,3)
 299       enddo
 300       write (2,'(a)') 'TER'
 301       do i=1,nnsup_t
 302         ii=isel(i)
 303         call mvvad(b,xyz2(1,ii),t,z(1,ii))
 304         write (2,100) ii+resnum3-1,atnam(ii),waby(ii),ires(ii),
 305      &                (z(j,ii),j=1,3)
 306       enddo
 307       write (2,'(a)') 'TER'
 308
 309       ind=0
 310       do i=1,nsup
 311         insel=.false.
 312         do j=1,nsup_t
 313           if (issel(i).eq.isel(j)) then
 314             insel=.true.
 315             exit
 316           endif
 317         enddo
 318         if (.not.insel) then
 319           ind=ind+1
 320           issel(ind)=issel(i)
 321         endif
 322       enddo
 323       nsup=nsup-nsup_t
 324       print *,"nsup",nsup
 325       if (nsup.gt.0) print *,"issel",(issel(i),i=1,nsup)
 326       do i=1,nsup
 327         do j=1,3
 328           s1(j,i)=ca1(j,issel(i))
 329           s2(j,i)=ca2(j,issel(i))
 330         enddo
 331       enddo
 332
 333       ENDDO
 334
 335       write (2,'(5000i2)') (frag(i),i=1,resnum2-resnum1+1)
 336       stop 777
 337   111 print '(2a)','Error opening input file ',fnam1(:len_trim(fnam1))
 338       stop 111
 339   112 print '(2a)','Error opening input file ',fnam2(:len_trim(fnam2))
 340       stop 112
 341   100 format ('ATOM',i7,2x,a2,2x,a3,i6,4x,3f8.3)
 342       end
 343 c----------------------------------------------------------------------
 344       subroutine distout(ncon,rms)
 345       parameter(ncol=10)
 346       dimension rms(200,200)
 347       dimension b(ncol)
 348       iout=2
 349 c      icant(i,j)=((ncon+ncon-j)*(j-1))/2+i-j
 350       write (iout,'(a)') 'The distance matrix'
 351       do 1 i=1,ncon,ncol
 352       nlim=min0(i+ncol-1,ncon)
 353       write (iout,1000) (k,k=i,nlim)
 354       write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
 355  1000 format (/8x,10(i4,3x))
 356  1020 format (/1x,80(1h-)/)
 357       do 2 j=i,ncon
 358       jlim=min0(j,nlim)
 359       if (jlim.eq.j) then
 360         b(jlim-i+1)=0.0d0
 361         jlim1=jlim-1
 362       else
 363         jlim1=jlim
 364       endif
 365       do 3 k=i,jlim1
 366     3 b(k-i+1)=rms(k,j)
 367       write (iout,1010) j,(b(k),k=1,jlim-i+1)
 368     2 continue
 369     1 continue
 370  1010 format (i5,3x,10(f6.2,1x))
 371       return
 372       end
 373 !--------------------------------------------------
 374       subroutine mysort(n, x, ipermut)
 375       implicit none
 376       real x(n),xsort(n)
 377       integer ipermut(n)
 378       integer i,j,imax,ipm,n
 379       real xtemp
 380       do i=1,n
 381         xtemp=x(i)
 382         imax=i
 383         do j=i+1,n
 384           if (x(j).lt.xtemp) then
 385             imax=j
 386             xtemp=x(j)
 387           endif
 388         enddo
 389         x(imax)=x(i)
 390         x(i)=xtemp
 391         ipm=ipermut(imax)
 392         ipermut(imax)=ipermut(i)
 393         ipermut(i)=ipm
 394       enddo
 395       return
 396       end
 397