source/maxlik/src_MD_T_maxlik-NEWCORR-PMF-5/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,ib,iprot,lprn)
   4       implicit none
   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       integer i,is,js,ishif1,jfrag,nc_match1_max,n_shif1,n_shif2,
  10      &  nc_req_set,istr,nc_match_max,nc_req,nc_match1,ib,iprot
  11       double precision nc_frac
  12       logical llocal,lprn
  13       nc_match_max=0
  14       do i=1,ncont_ref
  15         nc_match_max=nc_match_max+
  16      &   min0(icont_ref(2,i)-icont_ref(1,i)-1,3)
  17       enddo
  18       if (istr.eq.3) then
  19         nc_req=0
  20       else if (nc_req_set.eq.0) then
  21         nc_req=nc_match_max*nc_frac
  22       else
  23         nc_req = dmin1(nc_match_max*nc_frac+0.5d0,
  24      &    dfloat(nc_req_set)+1.0d-7)
  25       endif
  26 c      write (iout,*) "match_contact: nc_req:",nc_req
  27 c      write (iout,*) "nc_match_max",nc_match_max
  28 c      write (iout,*) "jfrag",jfrag," n_shif1",n_shif1,
  29 c     &   " n_shif2",n_shif2
  30 C Match current contact map against reference contact map; exit, if at least
  31 C half of the contacts match
  32       call ncont_match(nc_match,nc_match1,0,0,ncont_ref,icont_ref,
  33      &    ncont,icont,jfrag,llocal,ib,iprot,lprn)
  34       nc_match1_max=nc_match1
  35       if (lprn .and. nc_match.gt.0) write (iout,*)
  36      &  "Shift:",0,0," nc_match1",nc_match1,
  37      &  " nc_match=",nc_match," req'd",nc_req
  38       if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
  39      &    nc_req.eq.0 .and. nc_match.eq.1) then
  40          ishif1=0
  41          ishif2=0
  42          return
  43       endif
  44 C If sufficient matches are not found, try to shift contact maps up to three
  45 C positions.
  46       if (n_shif1.gt.0) then
  47       do is=1,n_shif1
  48 C The following four tries help to find shifted beta-sheet patterns
  49 C Shift "left" strand backward
  50         call ncont_match(nc_match,nc_match1,-is,0,ncont_ref,
  51      &    icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
  52         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
  53         if (lprn .and. nc_match.gt.0) write (iout,*)
  54      &    "Shift:",-is,0," nc_match1",nc_match1,
  55      &    " nc_match=",nc_match," req'd",nc_req
  56         if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
  57      &     nc_req.eq.0 .and. nc_match.eq.1) then
  58           ishif1=-is
  59           ishif2=0
  60           return
  61         endif
  62 C Shift "left" strand forward
  63         call ncont_match(nc_match,nc_match1,is,0,ncont_ref,
  64      &      icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
  65         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
  66         if (lprn .and. nc_match.gt.0) write (iout,*)
  67      &   "Shift:",is,0," nc_match1",nc_match1,
  68      &   " nc_match=",nc_match," req'd",nc_req
  69         if (nc_req.gt.0 .and. nc_match.ge.nc_req .or.
  70      &     nc_req.eq.0 .and. nc_match.eq.1) then
  71           ishif1=is
  72           ishif2=0
  73           return
  74         endif
  75       enddo
  76       if (nc_req.eq.0) return
  77 C Shift "right" strand backward
  78       do is=1,n_shif1
  79         call ncont_match(nc_match,nc_match1,0,-is,ncont_ref,
  80      &     icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
  81         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
  82         if (lprn .and. nc_match.gt.0) write (iout,*)
  83      &    "Shift:",0,-is," nc_match1",nc_match1,
  84      &    " nc_match=",nc_match," req'd",nc_req
  85         if (nc_match.ge.nc_req) then
  86           ishif1=0
  87           ishif2=-is
  88           return
  89         endif
  90 C Shift "right" strand upward
  91         call ncont_match(nc_match,nc_match1,0,is,ncont_ref,
  92      &    icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
  93         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
  94         if (lprn .and. nc_match.gt.0) write (iout,*)
  95      &    "Shift:",0,is," nc_match1",nc_match1,
  96      &    " nc_match=",nc_match," req'd",nc_req
  97         if (nc_match.ge.nc_req) then
  98           ishif1=0
  99           ishif2=is
 100           return
 101         endif
 102       enddo ! is
 103 C Now try to shift both residues in contacts.
 104       do is=1,n_shif1
 105         do js=1,is
 106           if (js.ne.is) then
 107             call ncont_match(nc_match,nc_match1,-is,-js,ncont_ref,
 108      &        icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 109             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 110             if (lprn .and. nc_match.gt.0) write (iout,*)
 111      &         "Shift:",-is,-js," nc_match1",nc_match1,
 112      &         " nc_match=",nc_match," req'd",nc_req
 113             if (nc_match.ge.nc_req) then
 114               ishif1=-is
 115               ishif2=-js
 116               return
 117             endif
 118             call ncont_match(nc_match,nc_match1,is,js,ncont_ref,
 119      &        icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 120             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 121             if (lprn .and. nc_match.gt.0) write (iout,*)
 122      &        "Shift:",is,js," nc_match1",nc_match1,
 123      &        " nc_match=",nc_match," req'd",nc_req
 124             if (nc_match.ge.nc_req) then
 125               ishif1=is
 126               ishif2=js
 127               return
 128             endif
 129 c
 130             call ncont_match(nc_match,nc_match1,-js,-is,ncont_ref,
 131      &        icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 132             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 133             if (lprn .and. nc_match.gt.0) write (iout,*)
 134      &        "Shift:",-js,-is," nc_match1",nc_match1,
 135      &        " nc_match=",nc_match," req'd",nc_req
 136             if (nc_match.ge.nc_req) then
 137               ishif1=-js
 138               ishif2=-is
 139               return
 140             endif
 141 c
 142             call ncont_match(nc_match,nc_match1,js,is,ncont_ref,
 143      &        icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 144             if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 145             if (lprn .and. nc_match.gt.0) write (iout,*)
 146      &        "Shift:",js,is," nc_match1",nc_match1,
 147      &        " nc_match=",nc_match," req'd",nc_req
 148             if (nc_match.ge.nc_req) then
 149               ishif1=js
 150               ishif2=is
 151               return
 152             endif
 153           endif
 154 c
 155           if (is+js.le.n_shif1) then
 156           call ncont_match(nc_match,nc_match1,-is,js,ncont_ref,
 157      &      icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 158           if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 159           if (lprn .and. nc_match.gt.0) write (iout,*)
 160      &     "Shift:",-is,js," nc_match1",nc_match1,
 161      &     " nc_match=",nc_match," req'd",nc_req
 162           if (nc_match.ge.nc_req) then
 163             ishif1=-is
 164             ishif2=js
 165             return
 166           endif
 167 c
 168           call ncont_match(nc_match,nc_match1,js,-is,ncont_ref,
 169      &      icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 170           if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 171           if (lprn .and. nc_match.gt.0) write (iout,*)
 172      &     "Shift:",js,-is," nc_match1",nc_match1,
 173      &     " nc_match=",nc_match," req'd",nc_req
 174           if (nc_match.ge.nc_req) then
 175             ishif1=js
 176             ishif2=-is
 177             return
 178           endif
 179           endif
 180 c
 181         enddo !js
 182       enddo !is
 183       endif
 184
 185       if (n_shif2.gt.0) then
 186       do is=1,n_shif2
 187         call ncont_match(nc_match,nc_match1,-is,-is,ncont_ref,
 188      &    icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 189         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 190         if (lprn .and. nc_match.gt.0) write (iout,*)
 191      &     "Shift:",-is,-is," nc_match1",nc_match1,
 192      &     " nc_match=",nc_match," req'd",nc_req
 193         if (nc_match.ge.nc_req) then
 194           ishif1=-is
 195           ishif2=-is
 196           return
 197         endif
 198         call ncont_match(nc_match,nc_match1,is,is,ncont_ref,
 199      &    icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 200         if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1
 201         if (lprn .and. nc_match.gt.0) write (iout,*)
 202      &    "Shift:",is,is," nc_match1",nc_match1,
 203      &    " nc_match=",nc_match," req'd",nc_req
 204         if (nc_match.ge.nc_req) then
 205           ishif1=is
 206           ishif2=is
 207           return
 208         endif
 209       enddo
 210       endif
 211 C If this point is reached, the contact maps are different.
 212       nc_match=0
 213       ishif1=0
 214       ishif2=0
 215       return
 216       end
 217 c-------------------------------------------------------------------------
 218       subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2,
 219      &   ncont_ref,icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
 220       implicit none
 221       include 'DIMENSIONS'
 222       include 'DIMENSIONS.ZSCOPT'
 223       include 'DIMENSIONS.COMPAR'
 224       include 'COMMON.IOUNITS'
 225       include 'COMMON.INTERACT'
 226       include 'COMMON.GEO'
 227       include 'COMMON.COMPAR'
 228       integer nc_match,nc_match1,ishif1,jfrag,i,j,ic1,ic2,ib,iprot
 229       double precision diffang,fract
 230       logical llocal,lprn
 231       integer ncont_ref,icont_ref(2,maxcont),ncont,icont(2,maxcont),
 232      &   icont_match(2,maxcont),ishift,ishif2,nang_pair,
 233      &   iang_pair(2,maxres)
 234 C Compare the contact map against the reference contact map; they're stored
 235 C in ICONT and ICONT_REF, respectively. The current contact map can be shifted.
 236       if (lprn) write (iout,'(80(1h*))')
 237       nc_match=0
 238       nc_match1=0
 239 c Check the local structure by comparing dihedral angles.
 240       if (llocal .and. ncont_ref.eq.0) then
 241 c If there are no contacts just compare the dihedral angles and exit.
 242         call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag,ib,iprot),
 243      &    diffang,fract,ib,iprot,lprn)
 244         if (lprn) write (iout,*) "diffang:",diffang*rad2deg,
 245      &   " ang_cut:",ang_cut(jfrag,ib,iprot)*rad2deg," fract",fract
 246         if (diffang.le.ang_cut(jfrag,ib,iprot) .and.
 247      &   fract.ge.frac_min(jfrag,ib,iprot))
 248      &  then
 249           nc_match=1
 250         else
 251           nc_match=0
 252         endif
 253         return
 254       endif
 255       nang_pair=0
 256       do i=1,ncont
 257         ic1=icont(1,i)+ishif1
 258         ic2=icont(2,i)+ishif2
 259 c        write (iout,*) "i",i," ic1",ic1," ic2",ic2
 260         if (ic1.lt.nnt .or. ic2.gt.nct) goto 10
 261         do j=1,ncont_ref
 262           if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then
 263             nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3)
 264             nc_match1=nc_match1+1
 265             icont_match(1,nc_match1)=ic1
 266             icont_match(2,nc_match1)=ic2
 267 c            call add_angpair(icont(1,i),icont_ref(1,j),
 268 c     &         nang_pair,iang_pair)
 269 c            call add_angpair(icont(2,i),icont_ref(2,j),
 270 c     &         nang_pair,iang_pair)
 271             if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i),
 272      &       " match",icont_ref(1,j),icont_ref(2,j),
 273      &       " shifts",ishif1,ishif2
 274             goto 10
 275           endif
 276         enddo
 277    10   continue
 278       enddo
 279       if (lprn) then
 280         write (iout,*) "nc_match",nc_match," nc_match1",nc_match1
 281         write (iout,*) "icont_match"
 282         do i=1,nc_match1
 283           write (iout,*) icont_match(1,i),icont_match(2,i)
 284         enddo
 285       endif
 286       if (llocal .and. nc_match.gt.0) then
 287         call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn,
 288      &    ang_cut1(jfrag,ib,iprot),diffang,fract,ib,iprot)
 289         if (lprn) write (iout,*) "diffang:",diffang*rad2deg,
 290      &   " ang_cut:",ang_cut(jfrag,ib,iprot)*rad2deg,
 291      &   " ang_cut1",ang_cut1(jfrag,ib,iprot)
 292         if (diffang.gt.ang_cut(jfrag,ib,iprot)
 293      &    .or. fract.lt.frac_min(jfrag,ib,iprot)) nc_match=0
 294       endif
 295 c      if (nc_match.gt.0) then
 296 c        diffang = angnorm1(nang_pair,iang_pair,lprn)
 297 c        if (diffang.gt.ang_cut(jfrag,ib,iprot)) nc_match=0
 298 c      endif
 299       if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,
 300      &   " diffang",rad2deg*diffang," nc_match",nc_match
 301       return
 302       end
 303 c------------------------------------------------------------------------------
 304       subroutine match_secondary(jfrag,isecstr,nsec_match,nsec_nomatch,
 305      &  avesig,ib,iprot,lprn)
 306 c This subroutine compares the secondary structure (isecstr) of fragment jfrag
 307 c conformation considered to that of the reference conformation.
 308 c Returns the number of equivalent residues (nsec_match).
 309       implicit none
 310       include 'DIMENSIONS'
 311       include 'DIMENSIONS.ZSCOPT'
 312       include 'DIMENSIONS.COMPAR'
 313       include 'COMMON.IOUNITS'
 314       include 'COMMON.CHAIN'
 315       include 'COMMON.PEPTCONT'
 316       include 'COMMON.COMPAR'
 317       logical lprn
 318       integer isecstr(maxres)
 319       integer i,j,j1,k,jc,jc1,jfrag,nsec_match,nsec_nomatch,ib,iprot,
 320      &  npart,nsig
 321       double precision sumsig,sumsigs1,sumsigs2,avesig,signature,sig,
 322      &  strand_signature
 323       external signature,strand_signature
 324       npart = npiece(jfrag,1,ib,iprot)
 325       nsec_match=0
 326       nsec_nomatch=0
 327       if (lprn) then
 328         write (iout,*) "match_secondary jfrag",jfrag," ifrag",
 329      &        (ifrag(1,i,jfrag,ib,iprot),ifrag(2,i,jfrag,ib,iprot),
 330      &        i=1,npart)
 331         write (iout,'(80i1)') (isec_ref(j,iprot),j=1,nres)
 332         write (iout,'(80i1)') (isecstr(j),j=1,nres)
 333       endif
 334       do i=1,npart
 335         do j=ifrag(1,i,jfrag,ib,iprot),ifrag(2,i,jfrag,ib,iprot)
 336 c The residue has equivalent conformational state to that of the reference
 337 c structure, if:
 338 c  a) the conformational states are equal or
 339 c  b) the reference state is a coil and that of the conformation considered
 340 c     is a strand or
 341 c  c) the conformational state of the conformation considered is a strand
 342 c     and that of the reference conformation is a coil.
 343 c 10/28/02 - case (b) deleted.
 344           if (isecstr(j).eq.isec_ref(j,iprot) .or.
 345 c     &        isecstr(j).eq.0 .and. isec_ref(j,iprot).eq.1 .or.
 346      &        isec_ref(j,iprot).eq.0 .and. isecstr(j).eq.1)
 347      &      nsec_match=nsec_match+1
 348 c 8/12/05 Unmatched secondary structure
 349           if (isecstr(j).eq.1 .and. isec_ref(j,iprot).eq.2 .or.
 350      &        isecstr(j).eq.2 .and.
 351      &        (isec_ref(j,iprot).eq.0 .or. isec_ref(j,iprot).eq.1)
 352      &        .or. isecstr(j).eq.3 .and. isec_ref(j,iprot).ne.3)
 353      &       nsec_nomatch = nsec_nomatch+1
 354         enddo
 355       enddo
 356       avesig=strand_signature(jfrag,ib,iprot,lprn)
 357       return
 358       end
 359 c------------------------------------------------------------------------------
 360       double precision function strand_signature(jfrag,ib,iprot,lprn)
 361 c This subroutine compares the secondary structure (isecstr) of fragment jfrag
 362 c conformation considered to that of the reference conformation.
 363 c Returns the number of equivalent residues (nsec_match).
 364       implicit none
 365       include 'DIMENSIONS'
 366       include 'DIMENSIONS.ZSCOPT'
 367       include 'DIMENSIONS.COMPAR'
 368       include 'COMMON.IOUNITS'
 369       include 'COMMON.CHAIN'
 370       include 'COMMON.PEPTCONT'
 371       include 'COMMON.COMPAR'
 372       logical lprn
 373       integer isecstr(maxres)
 374       integer i,j,j1,k,jc,jc1,jfrag,nsec_match,nsec_nomatch,ib,iprot,
 375      &  npart,nsig
 376       double precision sumsig,sumsigs1,sumsigs2,avesig,signature,sig
 377       external signature
 378 c 12/15/07 Compute signature for beta structure
 379       nsig=0
 380       sumsig=0.0d0
 381       sumsigs1=0.0d0
 382       sumsigs2=0.0d0
 383       npart = npiece(jfrag,1,ib,iprot)
 384       do i=1,npart
 385         do j=ifrag(1,i,jfrag,ib,iprot)+1,ifrag(2,i,jfrag,ib,iprot)-1
 386           if (isec_ref(j,iprot).eq.1) then
 387             do k=1,ncont_pept_near(j,ib,iprot)
 388               jc=icont_pept_near(k,j,ib,iprot)
 389               call find_dir(j,jc,j1,jc1,ib,iprot)
 390               sig=signature(j,jc,j1,jc1,.false.)
 391               sumsig=sumsig+sig*sig_ref(k,j,ib,iprot)
 392               sumsigs1=sumsigs1+sig_ref(k,j,ib,iprot)**2
 393               sumsigs2=sumsigs2+sig*sig
 394               nsig=nsig+1
 395               if (lprn) then
 396                 write (iout,*) "j",j," jc",jc," j1",j1," jc1",jc1,
 397      &           " isec_ref",isec_ref(j,iprot)," sig",sig,
 398      &           " sig_ref",sig_ref(k,j,ib,iprot)," sumsig",sumsig,
 399      &           " nsig",nsig
 400               endif
 401             enddo
 402           endif
 403         enddo
 404       enddo
 405       if (nsig.gt.0) then
 406         avesig=sumsig/dsqrt(sumsigs1*sumsigs2)
 407       else
 408         avesig=1.0d0
 409       endif
 410       if (lprn) write (iout,*) "avesig",avesig
 411       strand_signature=avesig
 412       return
 413       end
 414 c------------------------------------------------------------------------------
 415       double precision function signature(j,jc,j1,jc1,lprn)
 416       implicit none
 417       include 'DIMENSIONS'
 418       include 'DIMENSIONS.ZSCOPT'
 419       include 'DIMENSIONS.COMPAR'
 420       include 'COMMON.IOUNITS'
 421       include 'COMMON.CHAIN'
 422       include 'COMMON.PEPTCONT'
 423       include 'COMMON.COMPAR'
 424       integer j,j1,jc,jc1,k
 425       double precision u(3),v(3),w(3),z(3),unorm,vnorm,znorm
 426       double precision scalar
 427       external scalar
 428       logical lprn
 429 c
 430 c This function computes the signature of a residue involved in a beta-sheet;
 431 c the signature is defined as the scalar product of a vector normal to the
 432 c beta-sheet surface at contact beteen residues j and jc and the unit vector
 433 c of the bisector of the CA(j-1)-CA(j)-CA(j+1) angle.
 434 c
 435 c A.Liwo 12/15/07
 436 c
 437 c Compute the unit vectors CA(j)-CA(jc) (u) and CA(j-1)-CA(j+1) (v)
 438 c
 439       do k=1,3
 440         u(k)=c(k,j1)+c(k,j+1)-c(k,jc1)-c(k,jc+1)
 441 c        u(k)=c(k,j)+c(k,j+1)-c(k,jc)-c(k,jc+1)
 442         v(k)=c(k,j+1)-c(k,j-1)
 443       enddo
 444       unorm=dsqrt(scalar(u,u))
 445       vnorm=dsqrt(scalar(v,v))
 446       do k=1,3
 447         u(k)=u(k)/unorm
 448         v(k)=v(k)/vnorm
 449       enddo
 450       if (lprn) then
 451         write (iout,*) "vector u",(u(k),k=1,3)
 452         write (iout,*) "vector v",(v(k),k=1,3)
 453       endif
 454 c
 455 c Compute the normal vector to the beta-sheet at residue j (defined as the
 456 c vector product of u and v).
 457 c
 458       call vecpr(u,v,w)
 459       if (lprn) write (iout,*) "vector w",(w(k),k=1,3)
 460 c
 461 c Compute the unit vector of the CA(j-1)-CA(j)-CA(j+1) bisector (z).
 462 c
 463       do k=1,3
 464         u(k)=c(k,j-1)-c(k,j)
 465         v(k)=c(k,j+1)-c(k,j)
 466       enddo
 467       unorm=dsqrt(scalar(u,u))
 468       vnorm=scalar(v,v)
 469       do k=1,3
 470         u(k)=u(k)/unorm
 471         v(k)=v(k)/vnorm
 472         z(k)=u(k)+v(k)
 473       enddo
 474       znorm=dsqrt(scalar(z,z))
 475       do k=1,3
 476         z(k)=z(k)/znorm
 477       enddo
 478       if (lprn) write (iout,*) "vector z",(z(k),k=1,3)
 479 c
 480 c Compute the signature
 481 c
 482       signature=scalar(w,z)
 483       return
 484       end
 485 c-------------------------------------------------------------------------------
 486       subroutine find_dir(j,jc,j1,jc1,ib,iprot)
 487       implicit none
 488       include 'DIMENSIONS'
 489       include 'DIMENSIONS.ZSCOPT'
 490       include 'DIMENSIONS.COMPAR'
 491       include 'COMMON.IOUNITS'
 492       include 'COMMON.CHAIN'
 493       include 'COMMON.PEPTCONT'
 494       include 'COMMON.COMPAR'
 495       integer j,j1,jc,jc1,l,ib,iprot
 496       double precision u(3),v(3),w(3),z(3),unorm,vnorm,znorm
 497       double precision scalar
 498       external scalar
 499       logical lprn
 500       j1=j
 501       jc1=jc
 502       do l=1,ncont_pept_near(j-1,ib,iprot)
 503         if (icont_pept_near(l,j-1,ib,iprot).eq.jc-1) then
 504           j1=j-1
 505           jc1=jc-1
 506           return
 507         else if (icont_pept_near(l,j-1,ib,iprot).eq.jc+1) then
 508           j1=j-1
 509           jc1=jc+2
 510           return
 511         endif
 512       enddo
 513       return
 514       end