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)
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
15 nc_match_max=nc_match_max+
16 & min0(icont_ref(2,i)-icont_ref(1,i)-1,3)
20 else if (nc_req_set.eq.0) then
21 nc_req=nc_match_max*nc_frac
23 nc_req = dmin1(nc_match_max*nc_frac+0.5d0,
24 & dfloat(nc_req_set)+1.0d-7)
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
44 C If sufficient matches are not found, try to shift contact maps up to three
46 if (n_shif1.gt.0) then
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
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
76 if (nc_req.eq.0) return
77 C Shift "right" strand backward
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
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
103 C Now try to shift both residues in contacts.
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
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
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
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
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
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
185 if (n_shif2.gt.0) then
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
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
211 C If this point is reached, the contact maps are different.
217 c-------------------------------------------------------------------------
218 subroutine ncont_match(nc_match,nc_match1,ishif1,ishif2,
219 & ncont_ref,icont_ref,ncont,icont,jfrag,llocal,ib,iprot,lprn)
222 include 'DIMENSIONS.ZSCOPT'
223 include 'DIMENSIONS.COMPAR'
224 include 'COMMON.IOUNITS'
225 include 'COMMON.INTERACT'
227 include 'COMMON.COMPAR'
228 integer nc_match,nc_match1,ishif1,jfrag,i,j,ic1,ic2,ib,iprot
229 double precision diffang,fract
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*))')
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))
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
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
280 write (iout,*) "nc_match",nc_match," nc_match1",nc_match1
281 write (iout,*) "icont_match"
283 write (iout,*) icont_match(1,i),icont_match(2,i)
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
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
299 if (lprn) write (iout,*) "ishif1",ishif1," ishif2",ishif2,
300 & " diffang",rad2deg*diffang," nc_match",nc_match
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).
311 include 'DIMENSIONS.ZSCOPT'
312 include 'DIMENSIONS.COMPAR'
313 include 'COMMON.IOUNITS'
314 include 'COMMON.CHAIN'
315 include 'COMMON.PEPTCONT'
316 include 'COMMON.COMPAR'
318 integer isecstr(maxres)
319 integer i,j,j1,k,jc,jc1,jfrag,nsec_match,nsec_nomatch,ib,iprot,
321 double precision sumsig,sumsigs1,sumsigs2,avesig,signature,sig,
323 external signature,strand_signature
324 npart = npiece(jfrag,1,ib,iprot)
328 write (iout,*) "match_secondary jfrag",jfrag," ifrag",
329 & (ifrag(1,i,jfrag,ib,iprot),ifrag(2,i,jfrag,ib,iprot),
331 write (iout,'(80i1)') (isec_ref(j,iprot),j=1,nres)
332 write (iout,'(80i1)') (isecstr(j),j=1,nres)
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
338 c a) the conformational states are equal or
339 c b) the reference state is a coil and that of the conformation considered
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
356 avesig=strand_signature(jfrag,ib,iprot,lprn)
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).
366 include 'DIMENSIONS.ZSCOPT'
367 include 'DIMENSIONS.COMPAR'
368 include 'COMMON.IOUNITS'
369 include 'COMMON.CHAIN'
370 include 'COMMON.PEPTCONT'
371 include 'COMMON.COMPAR'
373 integer isecstr(maxres)
374 integer i,j,j1,k,jc,jc1,jfrag,nsec_match,nsec_nomatch,ib,iprot,
376 double precision sumsig,sumsigs1,sumsigs2,avesig,signature,sig
378 c 12/15/07 Compute signature for beta structure
383 npart = npiece(jfrag,1,ib,iprot)
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
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,
406 avesig=sumsig/dsqrt(sumsigs1*sumsigs2)
410 if (lprn) write (iout,*) "avesig",avesig
411 strand_signature=avesig
414 c------------------------------------------------------------------------------
415 double precision function signature(j,jc,j1,jc1,lprn)
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
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.
437 c Compute the unit vectors CA(j)-CA(jc) (u) and CA(j-1)-CA(j+1) (v)
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)
444 unorm=dsqrt(scalar(u,u))
445 vnorm=dsqrt(scalar(v,v))
451 write (iout,*) "vector u",(u(k),k=1,3)
452 write (iout,*) "vector v",(v(k),k=1,3)
455 c Compute the normal vector to the beta-sheet at residue j (defined as the
456 c vector product of u and v).
459 if (lprn) write (iout,*) "vector w",(w(k),k=1,3)
461 c Compute the unit vector of the CA(j-1)-CA(j)-CA(j+1) bisector (z).
467 unorm=dsqrt(scalar(u,u))
474 znorm=dsqrt(scalar(z,z))
478 if (lprn) write (iout,*) "vector z",(z(k),k=1,3)
480 c Compute the signature
482 signature=scalar(w,z)
485 c-------------------------------------------------------------------------------
486 subroutine find_dir(j,jc,j1,jc1,ib,iprot)
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
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
507 else if (icont_pept_near(l,j-1,ib,iprot).eq.jc+1) then