--- /dev/null
+ subroutine contact(lprint,ncont,icont,co)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.FFIELD'
+ include 'COMMON.NAMES'
+ real*8 facont /1.569D0/ ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
+ integer ncont,icont(2,maxcont)
+ logical lprint
+ ncont=0
+ kkk=3
+ do i=nnt+kkk,nct
+ iti=iabs(itype(i))
+ do j=nnt,i-kkk
+ itj=iabs(itype(j))
+ if (ipot.ne.4) then
+c rcomp=sigmaii(iti,itj)+1.0D0
+ rcomp=facont*sigmaii(iti,itj)
+ else
+c rcomp=sigma(iti,itj)+1.0D0
+ rcomp=facont*sigma(iti,itj)
+ endif
+c rcomp=6.5D0
+c print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
+ if (dist(nres+i,nres+j).lt.rcomp) then
+ ncont=ncont+1
+ icont(1,ncont)=i
+ icont(2,ncont)=j
+ endif
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') 'Contact map:'
+ do i=1,ncont
+ i1=icont(1,i)
+ i2=icont(2,i)
+ it1=itype(i1)
+ it2=itype(i2)
+ write (iout,'(i3,2x,a,i4,2x,a,i4)')
+ & i,restyp(it1),i1,restyp(it2),i2
+ enddo
+ endif
+ co = 0.0d0
+ do i=1,ncont
+ co = co + dfloat(iabs(icont(1,i)-icont(2,i)))
+ enddo
+ co = co / (nres*ncont)
+ return
+ end
+c----------------------------------------------------------------------------
+ double precision function contact_fract(ncont,ncont_ref,
+ & icont,icont_ref)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
+ nmatch=0
+c print *,'ncont=',ncont,' ncont_ref=',ncont_ref
+c write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
+c write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
+c write (iout,'(20i4)') (icont(1,i),i=1,ncont)
+c write (iout,'(20i4)') (icont(2,i),i=1,ncont)
+ do i=1,ncont
+ do j=1,ncont_ref
+ if (icont(1,i).eq.icont_ref(1,j) .and.
+ & icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
+ enddo
+ enddo
+c print *,' nmatch=',nmatch
+c contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
+ contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
+ return
+ end
+c----------------------------------------------------------------------------
+ double precision function contact_fract_nn(ncont,ncont_ref,
+ & icont,icont_ref)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
+ nmatch=0
+c print *,'ncont=',ncont,' ncont_ref=',ncont_ref
+c write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
+c write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
+c write (iout,'(20i4)') (icont(1,i),i=1,ncont)
+c write (iout,'(20i4)') (icont(2,i),i=1,ncont)
+ do i=1,ncont
+ do j=1,ncont_ref
+ if (icont(1,i).eq.icont_ref(1,j) .and.
+ & icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
+ enddo
+ enddo
+c print *,' nmatch=',nmatch
+c contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
+ contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont)
+ return
+ end
+c----------------------------------------------------------------------------
+ subroutine hairpin(lprint,nharp,iharp)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.FFIELD'
+ include 'COMMON.NAMES'
+ integer ncont,icont(2,maxcont)
+ integer nharp,iharp(4,maxres/3)
+ logical lprint,not_done
+ real*8 rcomp /6.0d0/
+ ncont=0
+ kkk=0
+c print *,'nnt=',nnt,' nct=',nct
+ do i=nnt,nct-3
+ do k=1,3
+ c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
+ enddo
+ do j=i+2,nct-1
+ do k=1,3
+ c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
+ enddo
+ if (dist(2*nres+1,2*nres+2).lt.rcomp) then
+ ncont=ncont+1
+ icont(1,ncont)=i
+ icont(2,ncont)=j
+ endif
+ enddo
+ enddo
+ if (lprint) then
+ write (iout,'(a)') 'PP contact map:'
+ do i=1,ncont
+ i1=icont(1,i)
+ i2=icont(2,i)
+ it1=itype(i1)
+ it2=itype(i2)
+ write (iout,'(i3,2x,a,i4,2x,a,i4)')
+ & i,restyp(it1),i1,restyp(it2),i2
+ enddo
+ endif
+c finding hairpins
+ nharp=0
+ do i=1,ncont
+ i1=icont(1,i)
+ j1=icont(2,i)
+ if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then
+c write (iout,*) "found turn at ",i1,j1
+ ii1=i1
+ jj1=j1
+ not_done=.true.
+ do while (not_done)
+ i1=i1-1
+ j1=j1+1
+ do j=1,ncont
+ if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
+ enddo
+ not_done=.false.
+ 10 continue
+c write (iout,*) i1,j1,not_done
+ enddo
+ i1=i1+1
+ j1=j1-1
+ if (j1-i1.gt.4) then
+ nharp=nharp+1
+ iharp(1,nharp)=i1
+ iharp(2,nharp)=j1
+ iharp(3,nharp)=ii1
+ iharp(4,nharp)=jj1
+c write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4)
+ endif
+ endif
+ enddo
+c do i=1,nharp
+c write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
+c enddo
+ if (lprint) then
+ write (iout,*) "Hairpins:"
+ do i=1,nharp
+ i1=iharp(1,i)
+ j1=iharp(2,i)
+ ii1=iharp(3,i)
+ jj1=iharp(4,i)
+ write (iout,*)
+ write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=i1,ii1)
+ write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
+c do k=jj1,j1,-1
+c write (iout,'(a,i3,$)') restyp(itype(k)),k
+c enddo
+ enddo
+ endif
+ return
+ end
+c----------------------------------------------------------------------------
+
projects / unres.git / blobdiff