--- /dev/null
+ dimension s1(3,10000),s2(3,10000),ca1(3,10000),ca2(3,10000)
+ dimension xyz1(3,5000),xyz2(3,5000),waby(5000),ires(5000),
+ & atnam(5000),rms(500,500),imin(5000),imax(5000),isel(5000),
+ & isel_t(5000),issel(5000),diff(5000),iperm(5000)
+ character*4 atnam
+ dimension t(3),b(3,3),z(3,5000)
+ integer resnum,resnum1,resnum2,resnum3,resnum4,frag(5000)
+ character*80 karta
+ character*80 fnam1,fnam2,fnam3
+ logical insel
+ ifrag=1
+ nnn = iargc()
+ if (nnn.lt.3) then
+ stop
+ & 'Usage: suppdb pdb_file_1 pdb_file_2 out_file [resnum1 resnum2]'
+ endif
+ call getarg(1,fnam1)
+ call getarg(2,fnam2)
+ call getarg(3,fnam3)
+ open (2,file=fnam3,status='unknown')
+ if (nnn.lt.5) then
+ print *,'resnum1,resnum2'
+ read (*,*) resnum1,resnum2
+ else
+ call getarg(4,karta)
+ read (karta,*) resnum1
+ call getarg(5,karta)
+ read (karta,*) resnum2
+ endif
+ if (nnn.gt.5) then
+ call getarg(6,karta)
+ read (karta,*) resnum3
+ else
+ resnum3=resnum1
+ endif
+ resnum4=resnum3+resnum2-resnum1
+ print *,'resnum1',resnum1,' resnum2',resnum2
+ print *,'resnum3',resnum3,' resnum4',resnum4
+ print *,'opening file ',fnam1
+ open (1,file=fnam1,status='old',err=111)
+ iat=0
+ isup=0
+ write (2,'(/2a)') 'Coordinates of the first set from ',
+ & fnam1(:len_trim(fnam1))
+ do irec=1,10000
+ read (1,'(a80)',end=10) karta
+ if ((karta(1:4).eq.'ATOM'.or.karta(1:6).eq.'HETATM').and.
+ & karta(14:15).eq.'CA') then
+ read (karta(23:26),*) resnum
+ if (resnum.ge.resnum1 .and. resnum.le.resnum2) then
+ iat=iat+1
+ read (karta(14:16),'(a3)') atnam(iat)
+ read (karta(31:54),'(3f8.3)') xyz1(1,iat),xyz1(2,iat),
+ 1 xyz1(3,iat)
+ read (karta(23:26),'(i4)') ires(iat)
+ read (karta(18:20),'(a3)') waby(iat)
+ if (atnam(iat).eq.'CA ') then
+ write (2,'(a)') karta
+ isup=isup+1
+ do i=1,3
+ ca1(i,isup)=xyz1(i,iat)
+ s1(i,isup)=xyz1(i,iat)
+ enddo
+ endif
+ endif
+ endif
+ enddo
+ 10 nsup=isup
+ nat=iat
+ print *,'from first file: nat=',nat
+ close(1)
+ print *,'opening file ',fnam2
+ open (1,file=fnam2,status='old',err=112)
+ iat=0
+ isup=0
+ write (2,'(/2a)') 'Coordinates of the first set from ',
+ & fnam2(:len_trim(fnam2))
+ do irec=1,10000
+ read (1,'(a80)',end=20) karta
+ if ((karta(1:4).eq.'ATOM'.or.karta(1:6).eq.'HETATM').and.
+ & karta(14:15).eq.'CA') then
+ read (karta(23:26),*) resnum
+ if (resnum.ge.resnum3 .and. resnum.le.resnum4) then
+ iat=iat+1
+ read (karta(31:54),'(3f8.3)') xyz2(1,iat),xyz2(2,iat),
+ 1 xyz2(3,iat)
+ if (atnam(iat).eq.'CA ') then
+ write (2,'(a)') karta
+ isup=isup+1
+ do i=1,3
+ ca2(i,isup)=xyz2(i,iat)
+ s2(i,isup)=xyz2(i,iat)
+ enddo
+ endif
+ endif
+ endif
+ enddo
+ 20 continue
+ print *,'from second file: nat=',iat
+! write (2,'(//4a8)') '# length',' start',' rms'
+
+
+ do i=1,nsup
+ issel(i)=i
+ enddo
+
+ nsup_t = nsup
+
+ nsup_old=0
+
+ DO WHILE (NSUP>4 .and. nsup.ne.nsup_old)
+
+ nsup_old = nsup
+
+ do is1=4,nsup
+! write (2,'(a)')
+ do is2=1,nsup-is1+1
+ call fitsq(rms(is1,is2),s2(1,is2),s1(1,is2),is1,t,b)
+! write (2,'(2(i4,4x),2f8.2)') is1,is2+resnum1-1,rms(is1,is2)
+ enddo
+ enddo
+ write (2,'(//2a8,3a8)') '# length',' start','rmsmin',
+ & 'rmsmax','rmsave'
+ do i=4,nsup
+ imin(i)=1
+ imax(i)=1
+ rmin=rms(i,1)
+ rmax=rms(i,1)
+ rmsave=rms(i,1)**2
+ do j=2,nsup-i+1
+ if ( rms(i,j).lt.rmin) then
+ imin(i)=j
+ rmin=rms(i,j)
+ else if ( rms(i,j).gt.rmax) then
+ imax(i)=j
+ rmax=rms(i,j)
+ endif
+ rmsave=rmsave+rms(i,j)*rms(i,j)
+ enddo
+ rmsave=sqrt(rmsave/(nsup-i+1))
+ write (2,'(2(i4,4x),3(f8.2))') i,imin(i)+resnum1-1,
+ & rms(i,imin(i)), ! rms(i,imin(i))/sqrt(float(i)),
+ & rms(i,imax(i)), ! rms(i,imax(i))/sqrt(float(i)),
+ & rmsave ! ,rmsave/sqrt(float(i))
+ enddo
+ rmscut=4.0d0
+ rmscut1=1.0d0*rmscut
+ do i4=nsup,4,-1
+ if (rms(i4,imin(i4)).le.rmscut) exit
+ enddo
+ call fitsq(rrms,s2(1,imin(i4)),s1(1,imin(i4)),i4,t,b)
+ nsup_t=nsup
+ do i=1,nsup
+ isel(i)=issel(i)
+ enddo
+ rrms=1000.0d0
+c do while (rrms.gt.rmscut)
+ do while (rrms.gt.rmscut)
+ nnsup_t=0
+ std=0.0d0
+ diffmax=-10000.0d0
+ do i=1,nsup_t
+ ii=isel(i)
+ call mvvad(b,ca2(1,ii),t,z(1,i))
+ diffi=0.0d0
+ do j=1,3
+ diffi=diffi+(z(j,i)-ca1(j,ii))**2
+ enddo
+ std=std+diffi
+ diffi=sqrt(diffi)
+! write (2,*) i,ii,diffi
+ diff(i)=diffi
+ if (diffi.gt.diffmax) diffmax=diffi
+ enddo
+ do i=1,nsup_t
+ iperm(i)=i
+ enddo
+ write(2,*) 'diff'
+ write(2,*) (diff(i),i=1,nsup_t)
+ call mysort(nsup_t,diff,iperm)
+ write(2,*) 'diff'
+ write(2,*) (diff(i),i=1,nsup_t)
+ std=sqrt(std/nsup_t)
+ write (2,*) 'std=',std
+ write (*,*) 'std=',std
+ sstd=dmin1(1.5d0*std,0.9999d0*diffmax)
+! do i=1,nsup_t
+! ii=isel(i)
+! if (diff(i).lt.sstd .or. diff(i).lt.rmscut) then
+! nnsup_t=nnsup_t+1
+! isel(nnsup_t)=ii
+! endif
+! enddo
+ sumdiff=0.0d0
+ do i=1,nsup_t
+ ii=isel(iperm(i))
+c sumdiff=sumdiff+diff(i)**2
+c if (dsqrt(sumdiff/dfloat(i)).lt.rmscut1) then
+ if (diff(i).lt.rmscut1) then
+ nnsup_t=nnsup_t+1
+ isel_t(nnsup_t)=ii
+ endif
+ enddo
+ nsup_t=nnsup_t
+ do i=1,nsup_t
+ iperm(i)=i
+ isel(i)=isel_t(i)
+ enddo
+ if (nsup_t.eq.0 .or. nnsup_t.eq.0) exit
+ write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
+ write (2,'(20i4)') (isel(i),i=1,nnsup_t)
+ call mysort(nsup_t,isel,iperm)
+ print *,'nsup_t',nsup_t,' nnsup_t',nnsup_t
+ print '(20i4)',(isel(i),i=1,nnsup_t)
+ write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
+ write (2,'(20i4)') (isel(i),i=1,nnsup_t)
+ len_min=4
+ i=1
+ len=1
+ do i=1,nsup_t
+ ii=isel(i)
+ do j=1,3
+ s1(j,i)=ca1(j,ii)
+ s2(j,i)=ca2(j,ii)
+ enddo
+ enddo
+ do i=1,nsup_t
+ ii=isel(i)
+ do j=1,3
+ s1(j,i)=ca1(j,ii)
+ s2(j,i)=ca2(j,ii)
+ enddo
+ enddo
+ if (nsup_t.le.2) exit
+ print *,"before fitsq: nsup_r",nsup_t
+ call fitsq(rrms,s2(1,1),s1(1,1),nsup_t,t,b)
+ write (2,*) 'nsup_t=',nsup_t,' rrms=',rrms
+ write (*,*) 'nsup_t=',nsup_t,' rrms=',rrms
+ enddo
+ i=1
+ do while (i.le.nnsup_t)
+ write (*,*) 'i=',i,' isel',isel(i),' isel1',isel(i+1)
+ if (i.eq.nnsup_t .or. isel(i+1)-isel(i).gt.1) then
+ if (len.lt.len_min) then
+ write (*,'(20i4)') (isel(k),k=1,nnsup_t)
+ write (*,*) 'len=',len,' len_min=',len_min
+ do j=i+1,nnsup_t
+ isel(j-len)=isel(j)
+ enddo
+ nnsup_t=nnsup_t-len
+! write (2,'(20i4)') (isel(k),k=1,nnsup_t)
+ i=i-len+1
+ len=1
+ else
+ len=1
+ i=i+1
+ endif
+ else
+ len=len+1
+ i=i+1
+ endif
+! write (2,*) 'i=',i,' len=',len
+ enddo
+ print *,'final nsup_t',nsup_t
+ print *,"final isel"
+ print '(20i4)',(isel(i),i=1,nnsup_t)
+ write (2,*) 'final nsup_t',nsup_t
+ write (2,*) "final isel"
+ write (2,'(20i4)') (isel(i),i=1,nnsup_t)
+ if (nnsup_t.gt.10) then
+ do i=1,nnsup_t
+ frag(isel(i))=ifrag
+ enddo
+ ifrag=ifrag+1
+ endif
+ print *,"----- end of iteration"
+ write (2,*) 'nsup_t',nsup_t,' nnsup_t',nnsup_t
+ nsup_t=nnsup_t
+ write (2,'(/4a)') 'REMARK superimposed coordinates from files ',
+ &fnam1(:len_trim(fnam1)),' and ',fnam2(:len_trim(fnam2))
+ write (2,'(a,f6.2,a)')
+ & 'REMARK Calpha rms deviation:',rrms,
+ & ' angstroms'
+ write (2,'(a,I5)') 'REMARK number of superposed residues',nsup_t
+ i1=isel(1)
+ do i=2,nnsup_t
+ if (isel(i)-isel(i-1).gt.1) then
+ write (2,'(a,i4,a,i4,a,f5.1,a)')
+ & 'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1
+ i1=isel(i)
+ endif
+ enddo
+ write (2,'(a,i4,a,i4,a,f5.1,a)')
+ & 'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1
+ do i=1,nsup_t
+ ii=isel(i)
+ write (2,100) ii+resnum1-1,atnam(ii),waby(ii),ires(ii),
+ & (xyz1(j,ii),j=1,3)
+ enddo
+ write (2,'(a)') 'TER'
+ do i=1,nnsup_t
+ ii=isel(i)
+ call mvvad(b,xyz2(1,ii),t,z(1,ii))
+ write (2,100) ii+resnum3-1,atnam(ii),waby(ii),ires(ii),
+ & (z(j,ii),j=1,3)
+ enddo
+ write (2,'(a)') 'TER'
+
+ ind=0
+ do i=1,nsup
+ insel=.false.
+ do j=1,nsup_t
+ if (issel(i).eq.isel(j)) then
+ insel=.true.
+ exit
+ endif
+ enddo
+ if (.not.insel) then
+ ind=ind+1
+ issel(ind)=issel(i)
+ endif
+ enddo
+ nsup=nsup-nsup_t
+ print *,"nsup",nsup
+ if (nsup.gt.0) print *,"issel",(issel(i),i=1,nsup)
+ do i=1,nsup
+ do j=1,3
+ s1(j,i)=ca1(j,issel(i))
+ s2(j,i)=ca2(j,issel(i))
+ enddo
+ enddo
+
+ ENDDO
+
+ write (2,'(5000i1)') (frag(i),i=1,resnum2-resnum1+1)
+ stop 777
+ 111 print '(2a)','Error opening input file ',fnam1(:len_trim(fnam1))
+ stop 111
+ 112 print '(2a)','Error opening input file ',fnam2(:len_trim(fnam2))
+ stop 112
+ 100 format ('ATOM',i7,2x,a2,2x,a3,i6,4x,3f8.3)
+ end
+c----------------------------------------------------------------------
+ subroutine distout(ncon,rms)
+ parameter(ncol=10)
+ dimension rms(200,200)
+ dimension b(ncol)
+ iout=2
+c icant(i,j)=((ncon+ncon-j)*(j-1))/2+i-j
+ write (iout,'(a)') 'The distance matrix'
+ do 1 i=1,ncon,ncol
+ nlim=min0(i+ncol-1,ncon)
+ write (iout,1000) (k,k=i,nlim)
+ write (iout,'(8h--------,10a)') ('-------',k=i,nlim)
+ 1000 format (/8x,10(i4,3x))
+ 1020 format (/1x,80(1h-)/)
+ do 2 j=i,ncon
+ jlim=min0(j,nlim)
+ if (jlim.eq.j) then
+ b(jlim-i+1)=0.0d0
+ jlim1=jlim-1
+ else
+ jlim1=jlim
+ endif
+ do 3 k=i,jlim1
+ 3 b(k-i+1)=rms(k,j)
+ write (iout,1010) j,(b(k),k=1,jlim-i+1)
+ 2 continue
+ 1 continue
+ 1010 format (i5,3x,10(f6.2,1x))
+ return
+ end
+!--------------------------------------------------
+ subroutine mysort(n, x, ipermut)
+ implicit none
+ real x(n),xsort(n)
+ integer ipermut(n)
+ integer i,j,imax,ipm,n
+ real xtemp
+ do i=1,n
+ xtemp=x(i)
+ imax=i
+ do j=i+1,n
+ if (x(j).lt.xtemp) then
+ imax=j
+ xtemp=x(j)
+ endif
+ enddo
+ x(imax)=x(i)
+ x(i)=xtemp
+ ipm=ipermut(imax)
+ ipermut(imax)=ipermut(i)
+ ipermut(i)=ipm
+ enddo
+ return
+ end
+
projects / unres.git / blobdiff