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)
6 dimension t(3),b(3,3),z(3,5000)
7 integer resnum,resnum1,resnum2,resnum3,resnum4,frag(5000)
9 character*80 fnam1,fnam2,fnam3
15 & 'Usage: suppdb pdb_file_1 pdb_file_2 out_file [resnum1 resnum2]'
20 open (2,file=fnam3,status='unknown')
22 print *,'resnum1,resnum2'
23 read (*,*) resnum1,resnum2
26 read (karta,*) resnum1
28 read (karta,*) resnum2
32 read (karta,*) resnum3
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)
43 write (2,'(/2a)') 'Coordinates of the first set from ',
44 & fnam1(:len_trim(fnam1))
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
52 read (karta(14:16),'(a3)') atnam(iat)
53 read (karta(31:54),'(3f8.3)') xyz1(1,iat),xyz1(2,iat),
55 read (karta(23:26),'(i4)') ires(iat)
56 read (karta(18:20),'(a3)') waby(iat)
57 if (atnam(iat).eq.'CA ') then
61 ca1(i,isup)=xyz1(i,iat)
62 s1(i,isup)=xyz1(i,iat)
70 print *,'from first file: nat=',nat
72 print *,'opening file ',fnam2
73 open (1,file=fnam2,status='old',err=112)
76 write (2,'(/2a)') 'Coordinates of the first set from ',
77 & fnam2(:len_trim(fnam2))
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
85 read (karta(31:54),'(3f8.3)') xyz2(1,iat),xyz2(2,iat),
87 if (atnam(iat).eq.'CA ') then
91 ca2(i,isup)=xyz2(i,iat)
92 s2(i,isup)=xyz2(i,iat)
99 print *,'from second file: nat=',iat
100 ! write (2,'(//4a8)') '# length',' start',' rms'
111 DO WHILE (NSUP>4 .and. nsup.ne.nsup_old)
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)
122 write (2,'(//2a8,3a8)') '# length',' start','rmsmin',
131 if ( rms(i,j).lt.rmin) then
134 else if ( rms(i,j).gt.rmax) then
138 rmsave=rmsave+rms(i,j)*rms(i,j)
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))
149 if (rms(i4,imin(i4)).le.rmscut) exit
151 call fitsq(rrms,s2(1,imin(i4)),s1(1,imin(i4)),i4,t,b)
157 c do while (rrms.gt.rmscut)
158 do while (rrms.gt.rmscut)
164 call mvvad(b,ca2(1,ii),t,z(1,i))
167 diffi=diffi+(z(j,i)-ca1(j,ii))**2
171 ! write (2,*) i,ii,diffi
173 if (diffi.gt.diffmax) diffmax=diffi
179 write(2,*) (diff(i),i=1,nsup_t)
180 call mysort(nsup_t,diff,iperm)
182 write(2,*) (diff(i),i=1,nsup_t)
184 write (2,*) 'std=',std
185 write (*,*) 'std=',std
186 sstd=dmin1(1.5d0*std,0.9999d0*diffmax)
189 ! if (diff(i).lt.sstd .or. diff(i).lt.rmscut) then
197 c sumdiff=sumdiff+diff(i)**2
198 c if (dsqrt(sumdiff/dfloat(i)).lt.rmscut1) then
199 if (diff(i).lt.rmscut1) then
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)
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
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
251 ! write (2,'(20i4)') (isel(k),k=1,nnsup_t)
262 ! write (2,*) 'i=',i,' len=',len
264 print *,'final nsup_t',nsup_t
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
276 print *,"----- end of iteration"
277 write (2,*) 'nsup_t',nsup_t,' nnsup_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,
284 write (2,'(a,I5)') 'REMARK number of superposed residues',nsup_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
293 write (2,'(a,i4,a,i4,a,f5.1,a)')
294 & 'REMARK fragment',i1+resnum1-1,' - ',isel(i-1)+resnum1-1
297 write (2,100) ii+resnum1-1,atnam(ii),waby(ii),ires(ii),
300 write (2,'(a)') 'TER'
303 call mvvad(b,xyz2(1,ii),t,z(1,ii))
304 write (2,100) ii+resnum3-1,atnam(ii),waby(ii),ires(ii),
307 write (2,'(a)') 'TER'
313 if (issel(i).eq.isel(j)) then
325 if (nsup.gt.0) print *,"issel",(issel(i),i=1,nsup)
328 s1(j,i)=ca1(j,issel(i))
329 s2(j,i)=ca2(j,issel(i))
335 write (2,'(5000i1)') (frag(i),i=1,resnum2-resnum1+1)
337 111 print '(2a)','Error opening input file ',fnam1(:len_trim(fnam1))
339 112 print '(2a)','Error opening input file ',fnam2(:len_trim(fnam2))
341 100 format ('ATOM',i7,2x,a2,2x,a3,i6,4x,3f8.3)
343 c----------------------------------------------------------------------
344 subroutine distout(ncon,rms)
346 dimension rms(200,200)
349 c icant(i,j)=((ncon+ncon-j)*(j-1))/2+i-j
350 write (iout,'(a)') 'The distance matrix'
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-)/)
367 write (iout,1010) j,(b(k),k=1,jlim-i+1)
370 1010 format (i5,3x,10(f6.2,1x))
373 !--------------------------------------------------
374 subroutine mysort(n, x, ipermut)
378 integer i,j,imax,ipm,n
384 if (x(j).lt.xtemp) then
392 ipermut(imax)=ipermut(i)