source/unres/src_MD-M/contact.f

   1       subroutine contact(lprint,ncont,icont,co)
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.IOUNITS'
   5       include 'COMMON.CHAIN'
   6       include 'COMMON.INTERACT'
   7       include 'COMMON.FFIELD'
   8       include 'COMMON.NAMES'
   9       real*8 facont /1.569D0/  ! facont = (2/(1-sqrt(1-1/4)))**(1/6)
  10       integer ncont,icont(2,maxcont)
  11       logical lprint
  12       ncont=0
  13       kkk=3
  14       do i=nnt+kkk,nct
  15         iti=itype(i)
  16         do j=nnt,i-kkk
  17           itj=itype(j)
  18           if (ipot.ne.4) then
  19 c           rcomp=sigmaii(iti,itj)+1.0D0
  20             rcomp=facont*sigmaii(iti,itj)
  21           else
  22 c           rcomp=sigma(iti,itj)+1.0D0
  23             rcomp=facont*sigma(iti,itj)
  24           endif
  25 c         rcomp=6.5D0
  26 c         print *,'rcomp=',rcomp,' dist=',dist(nres+i,nres+j)
  27           if (dist(nres+i,nres+j).lt.rcomp) then
  28             ncont=ncont+1
  29             icont(1,ncont)=i
  30             icont(2,ncont)=j
  31           endif
  32         enddo
  33       enddo
  34       if (lprint) then
  35         write (iout,'(a)') 'Contact map:'
  36         do i=1,ncont
  37           i1=icont(1,i)
  38           i2=icont(2,i)
  39           it1=itype(i1)
  40           it2=itype(i2)
  41           write (iout,'(i3,2x,a,i4,2x,a,i4)')
  42      &     i,restyp(it1),i1,restyp(it2),i2
  43         enddo
  44       endif
  45       co = 0.0d0
  46       do i=1,ncont
  47         co = co + dfloat(iabs(icont(1,i)-icont(2,i)))
  48       enddo
  49       co = co / (nres*ncont)
  50       return
  51       end
  52 c----------------------------------------------------------------------------
  53       double precision function contact_fract(ncont,ncont_ref,
  54      &                                     icont,icont_ref)
  55       implicit real*8 (a-h,o-z)
  56       include 'DIMENSIONS'
  57       include 'COMMON.IOUNITS'
  58       integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
  59       nmatch=0
  60 c     print *,'ncont=',ncont,' ncont_ref=',ncont_ref
  61 c     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
  62 c     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
  63 c     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
  64 c     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
  65       do i=1,ncont
  66         do j=1,ncont_ref
  67           if (icont(1,i).eq.icont_ref(1,j) .and.
  68      &        icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
  69         enddo
  70       enddo
  71 c     print *,' nmatch=',nmatch
  72 c     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
  73       contact_fract=dfloat(nmatch)/dfloat(ncont_ref)
  74       return
  75       end
  76 c----------------------------------------------------------------------------
  77       double precision function contact_fract_nn(ncont,ncont_ref,
  78      &                                     icont,icont_ref)
  79       implicit real*8 (a-h,o-z)
  80       include 'DIMENSIONS'
  81       include 'COMMON.IOUNITS'
  82       integer ncont,ncont_ref,icont(2,maxcont),icont_ref(2,maxcont)
  83       nmatch=0
  84 c     print *,'ncont=',ncont,' ncont_ref=',ncont_ref
  85 c     write (iout,'(20i4)') (icont_ref(1,i),i=1,ncont_ref)
  86 c     write (iout,'(20i4)') (icont_ref(2,i),i=1,ncont_ref)
  87 c     write (iout,'(20i4)') (icont(1,i),i=1,ncont)
  88 c     write (iout,'(20i4)') (icont(2,i),i=1,ncont)
  89       do i=1,ncont
  90         do j=1,ncont_ref
  91           if (icont(1,i).eq.icont_ref(1,j) .and.
  92      &        icont(2,i).eq.icont_ref(2,j)) nmatch=nmatch+1
  93         enddo
  94       enddo
  95 c     print *,' nmatch=',nmatch
  96 c     contact_fract=dfloat(nmatch)/dfloat(max0(ncont,ncont_ref))
  97       contact_fract_nn=dfloat(ncont-nmatch)/dfloat(ncont)
  98       return
  99       end
 100 c----------------------------------------------------------------------------
 101       subroutine hairpin(lprint,nharp,iharp)
 102       implicit real*8 (a-h,o-z)
 103       include 'DIMENSIONS'
 104       include 'COMMON.IOUNITS'
 105       include 'COMMON.CHAIN'
 106       include 'COMMON.INTERACT'
 107       include 'COMMON.FFIELD'
 108       include 'COMMON.NAMES'
 109       integer ncont,icont(2,maxcont)
 110       integer nharp,iharp(4,maxres/3)
 111       logical lprint,not_done
 112       real*8 rcomp /6.0d0/
 113       ncont=0
 114       kkk=0
 115 c     print *,'nnt=',nnt,' nct=',nct
 116       do i=nnt,nct-3
 117         do k=1,3
 118           c(k,2*nres+1)=0.5d0*(c(k,i)+c(k,i+1))
 119         enddo
 120         do j=i+2,nct-1
 121           do k=1,3
 122             c(k,2*nres+2)=0.5d0*(c(k,j)+c(k,j+1))
 123           enddo
 124           if (dist(2*nres+1,2*nres+2).lt.rcomp) then
 125             ncont=ncont+1
 126             icont(1,ncont)=i
 127             icont(2,ncont)=j
 128           endif
 129         enddo
 130       enddo
 131       if (lprint) then
 132         write (iout,'(a)') 'PP contact map:'
 133         do i=1,ncont
 134           i1=icont(1,i)
 135           i2=icont(2,i)
 136           it1=itype(i1)
 137           it2=itype(i2)
 138           write (iout,'(i3,2x,a,i4,2x,a,i4)')
 139      &     i,restyp(it1),i1,restyp(it2),i2
 140         enddo
 141       endif
 142 c finding hairpins
 143       nharp=0
 144       do i=1,ncont
 145         i1=icont(1,i)
 146         j1=icont(2,i)
 147         if (j1.eq.i1+2 .and. i1.gt.nnt .and. j1.lt.nct) then
 148 c          write (iout,*) "found turn at ",i1,j1
 149           ii1=i1
 150           jj1=j1
 151           not_done=.true.
 152           do while (not_done)
 153             i1=i1-1
 154             j1=j1+1
 155             do j=1,ncont
 156               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
 157             enddo
 158             not_done=.false.
 159   10        continue
 160 c            write (iout,*) i1,j1,not_done
 161           enddo
 162           i1=i1+1
 163           j1=j1-1
 164           if (j1-i1.gt.4) then
 165             nharp=nharp+1
 166             iharp(1,nharp)=i1
 167             iharp(2,nharp)=j1
 168             iharp(3,nharp)=ii1
 169             iharp(4,nharp)=jj1
 170 c            write (iout,*)'nharp',nharp,' iharp',(iharp(k,nharp),k=1,4)
 171           endif
 172         endif
 173       enddo
 174 c      do i=1,nharp
 175 c            write (iout,*)'i',i,' iharp',(iharp(k,i),k=1,4)
 176 c      enddo
 177       if (lprint) then
 178       write (iout,*) "Hairpins:"
 179       do i=1,nharp
 180         i1=iharp(1,i)
 181         j1=iharp(2,i)
 182         ii1=iharp(3,i)
 183         jj1=iharp(4,i)
 184         write (iout,*)
 185         write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=i1,ii1)
 186         write (iout,'(20(a,i3,1x))') (restyp(itype(k)),k,k=j1,jj1,-1)
 187 c        do k=jj1,j1,-1
 188 c         write (iout,'(a,i3,$)') restyp(itype(k)),k
 189 c        enddo
 190       enddo
 191       endif
 192       return
 193       end
 194 c----------------------------------------------------------------------------
 195