Fixed the homology gradient in finegrain mode
[unres.git] / source / wham / src / 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