X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Ffragments%2Fsuppdb_ca_dis.f;fp=source%2Ffragments%2Fsuppdb_ca_dis.f;h=50b09cef157af830b9d5b2b4070559abd93f10a6;hb=ee6fc3238c648029c75d83863fe62448d5af9358;hp=0000000000000000000000000000000000000000;hpb=bac661b05a651a029943830167e1984e78b8ca2d;p=unres.git

diff --git a/source/fragments/suppdb_ca_dis.f b/source/fragments/suppdb_ca_dis.f
new file mode 100644
index 0000000..50b09ce
--- /dev/null
+++ b/source/fragments/suppdb_ca_dis.f
@@ -0,0 +1,397 @@
+      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
+