X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC-NEWCORR%2Fmatch_contact.f;fp=source%2Fwham%2Fsrc-NEWSC-NEWCORR%2Fmatch_contact.f;h=3ec20368e527a4e850a208a03bb59fbfebc1c6b2;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git diff --git a/source/wham/src-NEWSC-NEWCORR/match_contact.f b/source/wham/src-NEWSC-NEWCORR/match_contact.f new file mode 100644 index 0000000..3ec2036 --- /dev/null +++ b/source/wham/src-NEWSC-NEWCORR/match_contact.f @@ -0,0 +1,339 @@ + subroutine match_contact(ishif1,ishif2,nc_match,nc_match1_max, + & ncont_ref,icont_ref,ncont,icont,jfrag,n_shif1,n_shif2, + & nc_frac,nc_req_set,istr,llocal,lprn) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'COMMON.IOUNITS' + integer ncont_ref,icont_ref(2,maxcont),ncont,icont(2,maxcont), + & ishift,ishif2,nc_match + double precision nc_frac + logical llocal,lprn + nc_match_max=0 + do i=1,ncont_ref + nc_match_max=nc_match_max+ + & min0(icont_ref(2,i)-icont_ref(1,i)-1,3) + enddo + if (istr.eq.3) then + nc_req=0 + else if (nc_req_set.eq.0) then + nc_req=nc_match_max*nc_frac + else + nc_req = dmin1(nc_match_max*nc_frac+0.5d0, + & dfloat(nc_req_set)+1.0d-7) + endif +c write (iout,*) "match_contact: nc_req:",nc_req +c write (iout,*) "nc_match_max",nc_match_max +c write (iout,*) "jfrag",jfrag," n_shif1",n_shif1, +c & " n_shif2",n_shif2 +C Match current contact map against reference contact map; exit, if at least +C half of the contacts match + call ncont_match(nc_match,nc_match1,0,0,ncont_ref,icont_ref, + & ncont,icont,jfrag,llocal,lprn) + nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",0,0," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. + & nc_req.eq.0 .and. nc_match.eq.1) then + ishif1=0 + ishif2=0 + return + endif +C If sufficient matches are not found, try to shift contact maps up to three +C positions. + if (n_shif1.gt.0) then + do is=1,n_shif1 +C The following four tries help to find shifted beta-sheet patterns +C Shift "left" strand backward + call ncont_match(nc_match,nc_match1,-is,0,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",-is,0," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. + & nc_req.eq.0 .and. nc_match.eq.1) then + ishif1=-is + ishif2=0 + return + endif +C Shift "left" strand forward + call ncont_match(nc_match,nc_match1,is,0,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",is,0," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_req.gt.0 .and. nc_match.ge.nc_req .or. + & nc_req.eq.0 .and. nc_match.eq.1) then + ishif1=is + ishif2=0 + return + endif + enddo + if (nc_req.eq.0) return +C Shift "right" strand backward + do is=1,n_shif1 + call ncont_match(nc_match,nc_match1,0,-is,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",0,-is," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=0 + ishif2=-is + return + endif +C Shift "right" strand upward + call ncont_match(nc_match,nc_match1,0,is,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",0,is," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=0 + ishif2=is + return + endif + enddo ! is +C Now try to shift both residues in contacts. + do is=1,n_shif1 + do js=1,is + if (js.ne.is) then + call ncont_match(nc_match,nc_match1,-is,-js,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",-is,-js," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=-is + ishif2=-js + return + endif + call ncont_match(nc_match,nc_match1,is,js,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",is,js," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=is + ishif2=js + return + endif +c + call ncont_match(nc_match,nc_match1,-js,-is,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",-js,-is," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=-js + ishif2=-is + return + endif +c + call ncont_match(nc_match,nc_match1,js,is,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",js,is," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=js + ishif2=is + return + endif + endif +c + if (is+js.le.n_shif1) then + call ncont_match(nc_match,nc_match1,-is,js,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",-is,js," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=-is + ishif2=js + return + endif +c + call ncont_match(nc_match,nc_match1,js,-is,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",js,-is," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=js + ishif2=-is + return + endif + endif +c + enddo !js + enddo !is + endif + + if (n_shif2.gt.0) then + do is=1,n_shif2 + call ncont_match(nc_match,nc_match1,-is,-is,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",-is,-is," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=-is + ishif2=-is + return + endif + call ncont_match(nc_match,nc_match1,is,is,ncont_ref, + & icont_ref,ncont,icont,jfrag,llocal,lprn) + if (nc_match1.gt.nc_match1_max) nc_match1_max=nc_match1 + if (lprn .and. nc_match.gt.0) write (iout,*) + & "Shift:",is,is," nc_match1",nc_match1, + & " nc_match=",nc_match," req'd",nc_req + if (nc_match.ge.nc_req) then + ishif1=is + ishif2=is + return + endif + enddo + endif +C If this point is reached, the contact maps are different. + nc_match=0 + ishif1=0 + ishif2=0 + return + end +c------------------------------------------------------------------------- + subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2, + & ncont_ref,icont_ref,ncont,icont,jfrag,llocal,lprn) + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.COMPAR' + include 'COMMON.IOUNITS' + include 'COMMON.INTERACT' + include 'COMMON.GEO' + include 'COMMON.COMPAR' + logical llocal,lprn + integer ncont_ref,icont_ref(2,maxcont),ncont,icont(2,maxcont), + & icont_match(2,maxcont),ishift,ishif2,nang_pair, + & iang_pair(2,maxres) +C Compare the contact map against the reference contact map; they're stored +C in ICONT and ICONT_REF, respectively. The current contact map can be shifted. + if (lprn) write (iout,'(80(1h*))') + nc_match=0 + nc_match1=0 +c Check the local structure by comparing dihedral angles. +c write (iout,*) "ncont_match: ncont_ref",ncont_ref," llocal",llocal + if (llocal .and. ncont_ref.eq.0) then +c If there are no contacts just compare the dihedral angles and exit. + call angnorm(jfrag,ishif1,ishif2,ang_cut1(jfrag),diffang,fract, + & lprn) + if (lprn) write (iout,*) "diffang:",diffang*rad2deg, + & " ang_cut:",ang_cut(jfrag)*rad2deg," fract",fract + if (diffang.le.ang_cut(jfrag) .and. fract.ge.frac_min(jfrag)) + & then + nc_match=1 + else + nc_match=0 + endif + return + endif + nang_pair=0 + do i=1,ncont + ic1=icont(1,i)+ishif1 + ic2=icont(2,i)+ishif2 +c write (iout,*) "i",i," ic1",ic1," ic2",ic2 + if (ic1.lt.nnt .or. ic2.gt.nct) goto 10 + do j=1,ncont_ref + if (ic1.eq.icont_ref(1,j).and.ic2.eq.icont_ref(2,j)) then + nc_match=nc_match+min0(icont_ref(2,j)-icont_ref(1,j)-1,3) + nc_match1=nc_match1+1 + icont_match(1,nc_match1)=ic1 + icont_match(2,nc_match1)=ic2 +c call add_angpair(icont(1,i),icont_ref(1,j), +c & nang_pair,iang_pair) +c call add_angpair(icont(2,i),icont_ref(2,j), +c & nang_pair,iang_pair) + if (lprn) write (iout,*) "Contacts:",icont(1,i),icont(2,i), + & " match",icont_ref(1,j),icont_ref(2,j), + & " shifts",ishif1,ishif2 + goto 10 + endif + enddo + 10 continue + enddo + if (lprn) then + write (iout,*) "nc_match",nc_match," nc_match1",nc_match1 + write (iout,*) "icont_match" + do i=1,nc_match1 + write (iout,*) icont_match(1,i),icont_match(2,i) + enddo + endif + if (llocal .and. nc_match.gt.0) then + call angnorm2(jfrag,ishif1,ishif2,nc_match1,icont_match,lprn, + & ang_cut1(jfrag),diffang,fract) + if (lprn) write (iout,*) "diffang:",diffang*rad2deg, + & " ang_cut:",ang_cut(jfrag)*rad2deg, + & " ang_cut1",ang_cut1(jfrag)*rad2deg + if (diffang.gt.ang_cut(jfrag) + & .or. fract.lt.frac_min(jfrag)) nc_match=0 + endif +c if (nc_match.gt.0) then +c diffang = angnorm1(nang_pair,iang_pair,lprn) +c if (diffang.gt.ang_cut(jfrag)) nc_match=0 +c endif + if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2, + & " diffang",rad2deg*diffang," nc_match",nc_match + return + end +c------------------------------------------------------------------------------ + subroutine match_secondary(jfrag,isecstr,nsec_match,lprn) +c This subroutine compares the secondary structure (isecstr) of fragment jfrag +c conformation considered to that of the reference conformation. +c Returns the number of equivalent residues (nsec_match). + implicit real*8 (a-h,o-z) + include 'DIMENSIONS' + include 'DIMENSIONS.ZSCOPT' + include 'DIMENSIONS.COMPAR' + include 'COMMON.IOUNITS' + include 'COMMON.CHAIN' + include 'COMMON.PEPTCONT' + include 'COMMON.COMPAR' + logical lprn + integer isecstr(maxres) + npart = npiece(jfrag,1) + nsec_match=0 + if (lprn) then + write (iout,*) "match_secondary jfrag",jfrag," ifrag", + & (ifrag(1,i,jfrag),ifrag(2,i,jfrag),i=1,npart) + write (iout,'(80i1)') (isec_ref(j),j=1,nres) + write (iout,'(80i1)') (isecstr(j),j=1,nres) + endif + do i=1,npart + do j=ifrag(1,i,jfrag),ifrag(2,i,jfrag) +c The residue has equivalent conformational state to that of the reference +c structure, if: +c a) the conformational states are equal or +c b) the reference state is a coil and that of the conformation considered +c is a strand or +c c) the conformational state of the conformation considered is a strand +c and that of the reference conformation is a coil. +c 10/28/02 - case (b) deleted. + if (isecstr(j).eq.isec_ref(j) .or. +c & isecstr(j).eq.0 .and. isec_ref(j).eq.1 .or. + & isec_ref(j).eq.0 .and. isecstr(j).eq.1) + & nsec_match=nsec_match+1 + enddo + enddo + return + end