source/wham/src-M-Pawel/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,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
   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,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,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,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,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,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,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,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,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,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,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,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,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,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,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)
 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) write (iout,'(80(1h*))')
 233       nc_match=0
 234       nc_match1=0
 235 c Check the local structure by comparing dihedral angles.
 236 c      write (iout,*) "ncont_match: ncont_ref",ncont_ref," llocal",llocal
 237       if (llocal .and. ncont_ref.eq.0) then
 238 c If there are no contacts just compare the dihedral angles and exit.
 239         call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag),diffang,fract,
 240      &    lprn)
 241         if (lprn) write (iout,*) "diffang:",diffang*rad2deg,
 242      &   " ang_cut:",ang_cut(jfrag)*rad2deg," fract",fract
 243         if (diffang.le.ang_cut(jfrag) .and. fract.ge.frac_min(jfrag))
 244      &  then
 245           nc_match=1
 246         else
 247           nc_match=0
 248         endif
 249         return
 250       endif
 251       nang_pair=0
 252       do i=1,ncont
 253         ic1=icont(1,i)+ishif1
 254         ic2=icont(2,i)+ishif2
 255 c        write (iout,*) "i",i," ic1",ic1," ic2",ic2
 256         if (ic1.lt.nnt .or. ic2.gt.nct) goto 10
 257         do j=1,ncont_ref
 258           if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then
 259             nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3)
 260             nc_match1=nc_match1+1
 261             icont_match(1,nc_match1)=ic1
 262             icont_match(2,nc_match1)=ic2
 263 c            call add_angpair(icont(1,i),icont_ref(1,j),
 264 c     &         nang_pair,iang_pair)
 265 c            call add_angpair(icont(2,i),icont_ref(2,j),
 266 c     &         nang_pair,iang_pair)
 267             if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i),
 268      &       " match",icont_ref(1,j),icont_ref(2,j),
 269      &       " shifts",ishif1,ishif2
 270             goto 10
 271           endif
 272         enddo
 273    10   continue
 274       enddo
 275       if (lprn) then
 276         write (iout,*) "nc_match",nc_match," nc_match1",nc_match1
 277         write (iout,*) "icont_match"
 278         do i=1,nc_match1
 279           write (iout,*) icont_match(1,i),icont_match(2,i)
 280         enddo
 281       endif
 282       if (llocal .and. nc_match.gt.0) then
 283         call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn,
 284      &    ang_cut1(jfrag),diffang,fract)
 285         if (lprn) write (iout,*) "diffang:",diffang*rad2deg,
 286      &   " ang_cut:",ang_cut(jfrag)*rad2deg,
 287      &   " ang_cut1",ang_cut1(jfrag)*rad2deg
 288         if (diffang.gt.ang_cut(jfrag)
 289      &    .or. fract.lt.frac_min(jfrag)) nc_match=0
 290       endif
 291 c      if (nc_match.gt.0) then
 292 c        diffang = angnorm1(nang_pair,iang_pair,lprn)
 293 c        if (diffang.gt.ang_cut(jfrag)) nc_match=0
 294 c      endif
 295       if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,
 296      &   " diffang",rad2deg*diffang," nc_match",nc_match
 297       return
 298       end
 299 c------------------------------------------------------------------------------
 300       subroutine match_secondary(jfrag,isecstr,nsec_match,lprn)
 301 c This subroutine compares the secondary structure (isecstr) of fragment jfrag
 302 c conformation considered to that of the reference conformation.
 303 c Returns the number of equivalent residues (nsec_match).
 304       implicit real*8 (a-h,o-z)
 305       include 'DIMENSIONS'
 306       include 'DIMENSIONS.ZSCOPT'
 307       include 'DIMENSIONS.COMPAR'
 308       include 'COMMON.IOUNITS'
 309       include 'COMMON.CHAIN'
 310       include 'COMMON.PEPTCONT'
 311       include 'COMMON.COMPAR'
 312       logical lprn
 313       integer isecstr(maxres)
 314       npart = npiece(jfrag,1)
 315       nsec_match=0
 316       if (lprn) then
 317         write (iout,*) "match_secondary jfrag",jfrag," ifrag",
 318      &        (ifrag(1,i,jfrag),ifrag(2,i,jfrag),i=1,npart)
 319         write (iout,'(80i1)') (isec_ref(j),j=1,nres)
 320         write (iout,'(80i1)') (isecstr(j),j=1,nres)
 321       endif
 322       do i=1,npart
 323         do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag)
 324 c The residue has equivalent conformational state to that of the reference
 325 c structure, if:
 326 c  a) the conformational states are equal or
 327 c  b) the reference state is a coil and that of the conformation considered
 328 c     is a strand or
 329 c  c) the conformational state of the conformation considered is a strand
 330 c     and that of the reference conformation is a coil.
 331 c 10/28/02 - case (b) deleted.
 332           if (isecstr(j).eq.isec_ref(j) .or.
 333 c     &        isecstr(j).eq.0 .and. isec_ref(j).eq.1 .or.
 334      &        isec_ref(j).eq.0 .and. isecstr(j).eq.1)
 335      &      nsec_match=nsec_match+1
 336         enddo
 337       enddo
 338       return
 339       end