source/wham/src-HCD-5D/match_contact.f

   1       subroutine match_contact(ishif1,ishif2,nc_match,nc_match1_max,
   2      &   ncont_ref,icont_ref,ncont,icont,jfrag,n_shif1,n_shif2,
   3      &   nc_frac,nc_req_set,istr,ipermmin,llocal,lprn)
   4       implicit real*8 (a-h,o-z)
   5       include 'DIMENSIONS'
   6       include 'COMMON.IOUNITS'
   7       integer ncont_ref,icont_ref(2,maxcont),ncont,icont(2,maxcont),
   8      &   ishift,ishif2,nc_match,ipermmin
   9       double precision nc_frac
  10       logical llocal,lprn
  11       nc_match_max=0
  12       do i=1,ncont_ref
  13         nc_match_max=nc_match_max+
  14      &   min0(icont_ref(2,i)-icont_ref(1,i)-1,3)
  15       enddo
  16       if (istr.eq.3) then
  17         nc_req=0
  18       else if (nc_req_set.eq.0) then
  19         nc_req=nc_match_max*nc_frac
  20       else
  21         nc_req = dmin1(nc_match_max*nc_frac+0.5d0,
  22      &    dfloat(nc_req_set)+1.0d-7)
  23       endif
  24 c      write (iout,*) "match_contact: nc_req:",nc_req
  25 c      write (iout,*) "nc_match_max",nc_match_max
  26 c      write (iout,*) "jfrag",jfrag," n_shif1",n_shif1,
  27 c     &   " n_shif2",n_shif2
  28 C Match current contact map against reference contact map; exit, if at least
  29 C half of the contacts match
  30       call ncont_match(nc_match,nc_match1,0,0,ncont_ref,icont_ref,
  31      &    ncont,icont,jfrag,ipermmin,llocal,lprn)
  32       nc_match1_max=nc_match1
  33       if (lprn .and. nc_match.gt.0) write (iout,*)
  34      &  "Shift:",0,0," nc_match1",nc_match1,
  35      &  " nc_match=",nc_match," req'd",nc_req
  36       if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
  37      &    nc_req.eq.0 .and. nc_match.eq.1) then
  38          ishif1=0
  39          ishif2=0
  40          return
  41       endif
  42 C If sufficient matches are not found, try to shift contact maps up to three
  43 C positions.
  44       if (n_shif1.gt.0) then
  45       do is=1,n_shif1
  46 C The following four tries help to find shifted beta-sheet patterns
  47 C Shift "left" strand backward
  48         call ncont_match(nc_match,nc_match1,-is,0,ncont_ref,
  49      &    icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
  50         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
  51         if (lprn .and. nc_match.gt.0) write (iout,*)
  52      &    "Shift:",-is,0," nc_match1",nc_match1,
  53      &    " nc_match=",nc_match," req'd",nc_req
  54         if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
  55      &     nc_req.eq.0 .and. nc_match.eq.1) then
  56           ishif1=-is
  57           ishif2=0
  58           return
  59         endif
  60 C Shift "left" strand forward
  61         call ncont_match(nc_match,nc_match1,is,0,ncont_ref,
  62      &      icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
  63         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
  64         if (lprn .and. nc_match.gt.0) write (iout,*)
  65      &   "Shift:",is,0," nc_match1",nc_match1,
  66      &   " nc_match=",nc_match," req'd",nc_req
  67         if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
  68      &     nc_req.eq.0 .and. nc_match.eq.1) then
  69           ishif1=is
  70           ishif2=0
  71           return
  72         endif
  73       enddo
  74       if (nc_req.eq.0) return
  75 C Shift "right" strand backward
  76       do is=1,n_shif1
  77         call ncont_match(nc_match,nc_match1,0,-is,ncont_ref,
  78      &     icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
  79         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
  80         if (lprn .and. nc_match.gt.0) write (iout,*)
  81      &    "Shift:",0,-is," nc_match1",nc_match1,
  82      &    " nc_match=",nc_match," req'd",nc_req
  83         if (nc_match.ge.nc_req) then
  84           ishif1=0
  85           ishif2=-is
  86           return
  87         endif
  88 C Shift "right" strand upward
  89         call ncont_match(nc_match,nc_match1,0,is,ncont_ref,
  90      &    icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
  91         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
  92         if (lprn .and. nc_match.gt.0) write (iout,*)
  93      &    "Shift:",0,is," nc_match1",nc_match1,
  94      &    " nc_match=",nc_match," req'd",nc_req
  95         if (nc_match.ge.nc_req) then
  96           ishif1=0
  97           ishif2=is
  98           return
  99         endif
 100       enddo ! is
 101 C Now try to shift both residues in contacts.
 102       do is=1,n_shif1
 103         do js=1,is
 104           if (js.ne.is) then
 105             call ncont_match(nc_match,nc_match1,-is,-js,ncont_ref,
 106      &        icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 107             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 108             if (lprn .and. nc_match.gt.0) write (iout,*)
 109      &         "Shift:",-is,-js," nc_match1",nc_match1,
 110      &         " nc_match=",nc_match," req'd",nc_req
 111             if (nc_match.ge.nc_req) then
 112               ishif1=-is
 113               ishif2=-js
 114               return
 115             endif
 116             call ncont_match(nc_match,nc_match1,is,js,ncont_ref,
 117      &        icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 118             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 119             if (lprn .and. nc_match.gt.0) write (iout,*)
 120      &        "Shift:",is,js," nc_match1",nc_match1,
 121      &        " nc_match=",nc_match," req'd",nc_req
 122             if (nc_match.ge.nc_req) then
 123               ishif1=is
 124               ishif2=js
 125               return
 126             endif
 127 c
 128             call ncont_match(nc_match,nc_match1,-js,-is,ncont_ref,
 129      &        icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 130             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 131             if (lprn .and. nc_match.gt.0) write (iout,*)
 132      &        "Shift:",-js,-is," nc_match1",nc_match1,
 133      &        " nc_match=",nc_match," req'd",nc_req
 134             if (nc_match.ge.nc_req) then
 135               ishif1=-js
 136               ishif2=-is
 137               return
 138             endif
 139 c
 140             call ncont_match(nc_match,nc_match1,js,is,ncont_ref,
 141      &        icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 142             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 143             if (lprn .and. nc_match.gt.0) write (iout,*)
 144      &        "Shift:",js,is," nc_match1",nc_match1,
 145      &        " nc_match=",nc_match," req'd",nc_req
 146             if (nc_match.ge.nc_req) then
 147               ishif1=js
 148               ishif2=is
 149               return
 150             endif
 151           endif
 152 c
 153           if (is+js.le.n_shif1) then
 154           call ncont_match(nc_match,nc_match1,-is,js,ncont_ref,
 155      &      icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 156           if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 157           if (lprn .and. nc_match.gt.0) write (iout,*)
 158      &     "Shift:",-is,js," nc_match1",nc_match1,
 159      &     " nc_match=",nc_match," req'd",nc_req
 160           if (nc_match.ge.nc_req) then
 161             ishif1=-is
 162             ishif2=js
 163             return
 164           endif
 165 c
 166           call ncont_match(nc_match,nc_match1,js,-is,ncont_ref,
 167      &      icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 168           if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 169           if (lprn .and. nc_match.gt.0) write (iout,*)
 170      &     "Shift:",js,-is," nc_match1",nc_match1,
 171      &     " nc_match=",nc_match," req'd",nc_req
 172           if (nc_match.ge.nc_req) then
 173             ishif1=js
 174             ishif2=-is
 175             return
 176           endif
 177           endif
 178 c
 179         enddo !js
 180       enddo !is
 181       endif
 182
 183       if (n_shif2.gt.0) then
 184       do is=1,n_shif2
 185         call ncont_match(nc_match,nc_match1,-is,-is,ncont_ref,
 186      &    icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 187         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 188         if (lprn .and. nc_match.gt.0) write (iout,*)
 189      &     "Shift:",-is,-is," nc_match1",nc_match1,
 190      &     " nc_match=",nc_match," req'd",nc_req
 191         if (nc_match.ge.nc_req) then
 192           ishif1=-is
 193           ishif2=-is
 194           return
 195         endif
 196         call ncont_match(nc_match,nc_match1,is,is,ncont_ref,
 197      &    icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 198         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 199         if (lprn .and. nc_match.gt.0) write (iout,*)
 200      &    "Shift:",is,is," nc_match1",nc_match1,
 201      &    " nc_match=",nc_match," req'd",nc_req
 202         if (nc_match.ge.nc_req) then
 203           ishif1=is
 204           ishif2=is
 205           return
 206         endif
 207       enddo
 208       endif
 209 C If this point is reached, the contact maps are different.
 210       nc_match=0
 211       ishif1=0
 212       ishif2=0
 213       return
 214       end
 215 c-------------------------------------------------------------------------
 216       subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2,
 217      &   ncont_ref,icont_ref,ncont,icont,jfrag,ipermmin,llocal,lprn)
 218       implicit real*8 (a-h,o-z)
 219       include 'DIMENSIONS'
 220       include 'DIMENSIONS.ZSCOPT'
 221       include 'DIMENSIONS.COMPAR'
 222       include 'COMMON.IOUNITS'
 223       include 'COMMON.INTERACT'
 224       include 'COMMON.GEO'
 225       include 'COMMON.COMPAR'
 226       logical llocal,lprn
 227       integer ncont_ref,icont_ref(2,maxcont),ncont,icont(2,maxcont),
 228      &   icont_match(2,maxcont),ishift,ishif2,nang_pair,
 229      &   iang_pair(2,maxres),ipermmin,iperm
 230 C Compare the contact map against the reference contact map; they're stored
 231 C in ICONT and ICONT_REF, respectively. The current contact map can be shifted.
 232       if (lprn) then
 233          write (iout,'(80(1h*))')
 234          write (iout,*) "ncont_match"
 235          write (iout,*) "ipermmin",ipermmin
 236          write (iout,'(80(1h*))')
 237       endif
 238       nc_match=0
 239       nc_match1=0
 240 c Check the local structure by comparing dihedral angles.
 241 c      write (iout,*) "ncont_match: ncont_ref",ncont_ref," llocal",llocal
 242       if (llocal .and. ncont_ref.eq.0) then
 243 c If there are no contacts just compare the dihedral angles and exit.
 244         call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag),diffang,fract,
 245      &    ipermmin,lprn)
 246         if (lprn) write (iout,*) "diffang:",diffang*rad2deg,
 247      &   " ang_cut:",ang_cut(jfrag)*rad2deg," fract",fract
 248         if (diffang.le.ang_cut(jfrag) .and. fract.ge.frac_min(jfrag))
 249      &  then
 250           nc_match=1
 251         else
 252           nc_match=0
 253         endif
 254         return
 255       endif
 256       nang_pair=0
 257       do i=1,ncont
 258 c        write (iout,*) "i",i," icont",icont(1,i),icont(2,i)
 259         ic1=icont(1,i)+ishif1
 260         ic2=icont(2,i)+ishif2
 261 c        write (iout,*) "i",i," ic1",ic1," ic2",ic2
 262         if (ic1.lt.nnt .or. ic2.gt.nct) goto 10
 263         do j=1,ncont_ref
 264           if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then
 265             nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3)
 266             nc_match1=nc_match1+1
 267             icont_match(1,nc_match1)=ic1
 268             icont_match(2,nc_match1)=ic2
 269 c            call add_angpair(icont(1,i),icont_ref(1,j),
 270 c     &         nang_pair,iang_pair)
 271 c            call add_angpair(icont(2,i),icont_ref(2,j),
 272 c     &         nang_pair,iang_pair)
 273             if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i),
 274      &       " match",icont_ref(1,j),icont_ref(2,j),
 275      &       " shifts",ishif1,ishif2
 276             goto 10
 277           endif
 278         enddo
 279    10   continue
 280       enddo
 281       if (lprn) then
 282         write (iout,*) "nc_match",nc_match," nc_match1",nc_match1
 283         write (iout,*) "icont_match"
 284         do i=1,nc_match1
 285           write (iout,*) icont_match(1,i),icont_match(2,i)
 286         enddo
 287       endif
 288       if (llocal .and. nc_match.gt.0) then
 289         call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn,
 290      &    ang_cut1(jfrag),diffang,fract,ipermmin,lprn)
 291         if (lprn) write (iout,*) "diffang:",diffang*rad2deg,
 292      &   " ang_cut:",ang_cut(jfrag)*rad2deg,
 293      &   " ang_cut1",ang_cut1(jfrag)*rad2deg
 294         if (diffang.gt.ang_cut(jfrag)
 295      &    .or. fract.lt.frac_min(jfrag)) nc_match=0
 296       endif
 297 c      if (nc_match.gt.0) then
 298 c        diffang = angnorm1(nang_pair,iang_pair,lprn)
 299 c        if (diffang.gt.ang_cut(jfrag)) nc_match=0
 300 c      endif
 301       if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,
 302      &   " diffang",rad2deg*diffang," nc_match",nc_match
 303       return
 304       end
 305 c------------------------------------------------------------------------------
 306       subroutine match_secondary(jfrag,isecstr,nsec_match,ipermmin,lprn)
 307 c This subroutine compares the secondary structure (isecstr) of fragment jfrag
 308 c conformation considered to that of the reference conformation.
 309 c Returns the number of equivalent residues (nsec_match).
 310       implicit real*8 (a-h,o-z)
 311       include 'DIMENSIONS'
 312       include 'DIMENSIONS.ZSCOPT'
 313       include 'DIMENSIONS.COMPAR'
 314       include 'COMMON.IOUNITS'
 315       include 'COMMON.CHAIN'
 316       include 'COMMON.PEPTCONT'
 317       include 'COMMON.COMPAR'
 318       logical lprn
 319       integer isecstr(maxres),ipermmin,iperm
 320       npart = npiece(jfrag,1)
 321       nsec_match=0
 322       if (lprn) then
 323         write (iout,*) "match_secondary jfrag",jfrag," ifrag",
 324      &        (ifrag(1,i,jfrag),ifrag(2,i,jfrag),i=1,npart)
 325         write (iout,'(80i1)') (isec_ref(j),j=1,nres)
 326         write (iout,'(80i1)') (isecstr(j),j=1,nres)
 327       endif
 328       do i=1,npart
 329         do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
 330 c The residue has equivalent conformational state to that of the reference
 331 c structure, if:
 332 c  a) the conformational states are equal or
 333 c  b) the reference state is a coil and that of the conformation considered
 334 c     is a strand or
 335 c  c) the conformational state of the conformation considered is a strand
 336 c     and that of the reference conformation is a coil.
 337 c 10/28/02 - case (b) deleted.
 338           if (isecstr(iperm(j,ipermmin)).eq.isec_ref(j) .or.
 339 c     &        isecstr(j).eq.0 .and. isec_ref(j).eq.1 .or.
 340      &        isec_ref(j).eq.0 .and. isecstr(iperm(j,ipermmin)).eq.1)
 341      &      nsec_match=nsec_match+1
 342         enddo
 343       enddo
 344       return
 345       end