Merge branch 'devel' of mmka.chem.univ.gda.pl:unres into devel
[unres.git] / source / wham / src-M / secondary.f
1       subroutine define_fragments
2       implicit real*8 (a-h,o-z)
3       include 'DIMENSIONS'
4       include 'DIMENSIONS.ZSCOPT'
5       include 'DIMENSIONS.COMPAR'
6       include 'COMMON.IOUNITS'
7       include 'COMMON.TIME1'
8       include 'COMMON.FRAG'
9       include 'COMMON.SBRIDGE'
10       include 'COMMON.CONTROL'
11       include 'COMMON.COMPAR'
12       include 'COMMON.CHAIN'
13       include 'COMMON.HEADER'
14       include 'COMMON.GEO'
15       include 'COMMON.CONTACTS'
16       include 'COMMON.PEPTCONT'
17       include 'COMMON.INTERACT'
18       include 'COMMON.NAMES'
19       integer nstrand,istrand(2,maxres/2)
20       integer nhairp,ihairp(2,maxres/5) 
21       character*16 strstr(4) /'helix','hairpin','strand','strand pair'/
22       write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,
23      &               'NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,
24      &           'NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,
25      &  ' RMS_PAIR',irms_pair,' SPLIT_BET',isplit_bet
26       write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,
27      &  ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
28       write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,
29      &  ' MAXANG_HEL',angcut1_hel*rad2deg
30       write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,
31      &               ' MAXANG_BET',angcut1_bet*rad2deg
32       write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,
33      &               ' MAXANG_STRAND',angcut1_strand*rad2deg
34       write (iout,*) 'FRAC_MIN',frac_min_set
35 c Find secondary structure elements (helices and beta-sheets)
36       call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,
37      &   isec_ref)
38 c Define primary fragments. First include the helices.
39       nhairp=0
40       nstrand=0
41 c Merge helices
42 c AL 12/23/03 - to avoid splitting helices into very small fragments
43       if (merge_helices) then
44       write (iout,*) "Before merging helices: nhfrag",nhfrag
45       do i=1,nhfrag
46         write (2,*) hfrag(1,i),hfrag(2,i)
47       enddo
48       i=1
49       do while (i.lt.nhfrag)
50         if (hfrag(1,i+1)-hfrag(2,i).le.1) then
51           nhfrag=nhfrag-1
52           hfrag(2,i)=hfrag(2,i+1)
53           do j=i+1,nhfrag
54             hfrag(1,j)=hfrag(1,j+1)
55             hfrag(2,j)=hfrag(2,j+1)
56           enddo
57         endif 
58         i=i+1
59       enddo
60       write (iout,*) "After merging helices: nhfrag",nhfrag
61       do i=1,nhfrag
62         write (2,*) hfrag(1,i),hfrag(2,i)
63       enddo
64       endif
65       nfrag(1)=nhfrag
66       do i=1,nhfrag
67         npiece(i,1)=1
68         ifrag(1,1,i)=hfrag(1,i) 
69         ifrag(2,1,i)=hfrag(2,i) 
70         n_shift(1,i,1)=0
71         n_shift(2,i,1)=nshift_hel
72         ang_cut(i)=angcut_hel
73         ang_cut1(i)=angcut1_hel
74         frac_min(i)=frac_min_set
75         nc_fragm(i,1)=ncfrac_hel
76         nc_req_setf(i,1)=ncreq_hel
77         istruct(i)=1
78       enddo
79       write (iout,*) "isplit_bet",isplit_bet
80       if (isplit_bet.gt.1) then
81 c Split beta-sheets into strands and store strands as primary fragments.
82         call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
83         do i=1,nstrand
84           ii=i+nfrag(1)
85           npiece(ii,1)=1
86           ifrag(1,1,ii)=istrand(1,i)
87           ifrag(2,1,ii)=istrand(2,i)
88           n_shift(1,ii,1)=nshift_strand
89           n_shift(2,ii,1)=nshift_strand
90           ang_cut(ii)=angcut_strand
91           ang_cut1(ii)=angcut1_strand
92           frac_min(ii)=frac_min_set
93           nc_fragm(ii,1)=0
94           nc_req_setf(ii,1)=0
95           istruct(ii)=3
96         enddo
97         nfrag(1)=nfrag(1)+nstrand
98       else if (isplit_bet.eq.1) then
99 c Split only far beta-sheets; does not split hairpins.
100         call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
101         call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
102         do i=1,nhairp
103           ii=i+nfrag(1)
104           npiece(ii,1)=1
105           ifrag(1,1,ii)=ihairp(1,i) 
106           ifrag(2,1,ii)=ihairp(2,i) 
107           n_shift(1,ii,1)=nshift_bet
108           n_shift(2,ii,1)=nshift_bet
109           ang_cut(ii)=angcut_bet
110           ang_cut1(ii)=angcut1_bet
111           frac_min(ii)=frac_min_set
112           nc_fragm(ii,1)=ncfrac_bet
113           nc_req_setf(ii,1)=ncreq_bet
114           istruct(ii)=2
115         enddo
116         nfrag(1)=nfrag(1)+nhairp
117         do i=1,nstrand
118           ii=i+nfrag(1)
119           npiece(ii,1)=1
120           ifrag(1,1,ii)=istrand(1,i)
121           ifrag(2,1,ii)=istrand(2,i)
122           n_shift(1,ii,1)=nshift_strand
123           n_shift(2,ii,1)=nshift_strand
124           ang_cut(ii)=angcut_strand
125           ang_cut1(ii)=angcut1_strand
126           frac_min(ii)=frac_min_set
127           nc_fragm(ii,1)=0
128           nc_req_setf(ii,1)=0
129           istruct(ii)=3
130         enddo
131         nfrag(1)=nfrag(1)+nstrand
132       else
133 c Do not split beta-sheets; each pair of strands is a primary element.
134         call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
135         do i=1,nhairp
136           ii=i+nfrag(1)
137           npiece(ii,1)=1
138           ifrag(1,1,ii)=ihairp(1,i) 
139           ifrag(2,1,ii)=ihairp(2,i) 
140           n_shift(1,ii,1)=nshift_bet
141           n_shift(2,ii,1)=nshift_bet
142           ang_cut(ii)=angcut_bet
143           ang_cut1(ii)=angcut1_bet
144           frac_min(ii)=frac_min_set
145           nc_fragm(ii,1)=ncfrac_bet
146           nc_req_setf(ii,1)=ncreq_bet
147           istruct(ii)=2
148         enddo
149         nfrag(1)=nfrag(1)+nhairp
150         do i=1,nbfrag
151           ii=i+nfrag(1)
152           npiece(ii,1)=2
153           ifrag(1,1,ii)=bfrag(1,i) 
154           ifrag(2,1,ii)=bfrag(2,i) 
155           if (bfrag(3,i).lt.bfrag(4,i)) then
156             ifrag(1,2,ii)=bfrag(3,i)
157             ifrag(2,2,ii)=bfrag(4,i)
158           else
159             ifrag(1,2,ii)=bfrag(4,i)
160             ifrag(2,2,ii)=bfrag(3,i)
161           endif
162           n_shift(1,ii,1)=nshift_bet
163           n_shift(2,ii,1)=nshift_bet
164           ang_cut(ii)=angcut_bet
165           ang_cut1(ii)=angcut1_bet
166           frac_min(ii)=frac_min_set
167           nc_fragm(ii,1)=ncfrac_bet
168           nc_req_setf(ii,1)=ncreq_bet
169           istruct(ii)=4
170         enddo
171         nfrag(1)=nfrag(1)+nbfrag
172       endif
173       write (iout,*) "The following primary fragments were found:"
174       write (iout,*) "Helices:",nhfrag
175       do i=1,nhfrag
176         i1=ifrag(1,1,i)
177         i2=ifrag(2,1,i)
178         it1=itype(i1)
179         it2=itype(i2)
180         write (iout,'(i3,2x,a,i4,2x,a,i4)')
181      &       i,restyp(it1),i1,restyp(it2),i2
182       enddo
183       write (iout,*) "Hairpins:",nhairp
184       do i=nhfrag+1,nhfrag+nhairp
185         i1=ifrag(1,1,i)
186         i2=ifrag(2,1,i)
187         it1=itype(i1)
188         it2=itype(i2)
189         write (iout,'(i3,2x,a,i4,2x,a,i4,2x)')
190      &       i,restyp(it1),i1,restyp(it2),i2
191       enddo
192       write (iout,*) "Far strand pairs:",nbfrag
193       do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag
194         i1=ifrag(1,1,i)
195         i2=ifrag(2,1,i)
196         it1=itype(i1)
197         it2=itype(i2)
198         i3=ifrag(1,2,i)
199         i4=ifrag(2,2,i)
200         it3=itype(i3)
201         it4=itype(i4)
202         write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)')
203      &       i,restyp(it1),i1,restyp(it2),i2,
204      &         restyp(it3),i3,restyp(it4),i4
205       enddo
206       write (iout,*) "Strands:",nstrand
207       do i=nhfrag+nhairp+nbfrag+1,nfrag(1)
208         i1=ifrag(1,1,i)
209         i2=ifrag(2,1,i)
210         it1=itype(i1)
211         it2=itype(i2)
212         write (iout,'(i3,2x,a,i4,2x,a,i4)')
213      &       i,restyp(it1),i1,restyp(it2),i2
214       enddo
215       call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),
216      &  istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),
217      &  nc_fragm(1,1),nc_req_setf(1,1))
218       write (iout,*) "Fragments after sorting:"
219       do i=1,nfrag(1)
220         i1=ifrag(1,1,i)
221         i2=ifrag(2,1,i)
222         it1=itype(i1)
223         it2=itype(i2)
224         write (iout,'(i3,2x,a,i4,2x,a,i4,$)')
225      &       i,restyp(it1),i1,restyp(it2),i2
226         if (npiece(i,1).eq.1) then
227           write (iout,'(2x,a)') strstr(istruct(i))
228         else
229           i1=ifrag(1,2,i)
230           i2=ifrag(2,2,i)
231           it1=itype(i1)
232           it2=itype(i2)
233           write (iout,'(2x,a,i4,2x,a,i4,2x,a)')
234      &       restyp(it1),i1,restyp(it2),i2,strstr(istruct(i))
235         endif
236       enddo
237       return
238       end
239 c------------------------------------------------------------------------------
240       subroutine find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
241       implicit real*8 (a-h,o-z)
242       include 'DIMENSIONS'
243       include 'DIMENSIONS.ZSCOPT'
244       include 'DIMENSIONS.COMPAR'
245       include 'COMMON.IOUNITS'
246       integer nbfrag,bfrag(4,maxres/3)
247       integer nhairp,ihairp(2,maxres/5) 
248       write (iout,*) "Entered find_and_remove_hairpins"
249       write (iout,*) "nbfrag",nbfrag
250       do i=1,nbfrag
251         write (iout,*) i,(bfrag(k,i),k=1,4)
252       enddo
253       nhairp=0
254       i=1
255       do while (i.le.nbfrag)
256         write (iout,*) "check hairpin:",i,(bfrag(j,i),j=1,4)
257         if (bfrag(3,i).gt.bfrag(4,i) .and. bfrag(4,i)-bfrag(2,i).lt.5) 
258      &  then
259           write (iout,*) "Found hairpin:",i,bfrag(1,i),bfrag(3,i)
260           nhairp=nhairp+1
261           ihairp(1,nhairp)=bfrag(1,i)
262           ihairp(2,nhairp)=bfrag(3,i) 
263           nbfrag=nbfrag-1
264           do j=i,nbfrag
265             do k=1,4
266               bfrag(k,j)=bfrag(k,j+1)
267             enddo
268           enddo
269         else
270           i=i+1
271         endif
272       enddo
273       write (iout,*) "After finding hairpins:"
274       write (iout,*) "nhairp",nhairp
275       do i=1,nhairp
276         write (iout,*) i,ihairp(1,i),ihairp(2,i)
277       enddo
278       write (iout,*) "nbfrag",nbfrag
279       do i=1,nbfrag
280         write (iout,*) i,(bfrag(k,i),k=1,4)
281       enddo
282       return
283       end
284 c------------------------------------------------------------------------------
285       subroutine split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
286       implicit real*8 (a-h,o-z)
287       include 'DIMENSIONS'
288       include 'DIMENSIONS.ZSCOPT'
289       include 'DIMENSIONS.COMPAR'
290       include 'COMMON.IOUNITS'
291       integer nbfrag,bfrag(4,maxres/3)
292       integer nstrand,istrand(2,maxres/2)
293       integer nhairp,ihairp(2,maxres/5) 
294       logical found
295       write (iout,*) "Entered split_beta"
296       write (iout,*) "nbfrag",nbfrag
297       do i=1,nbfrag
298         write (iout,*) i,(bfrag(k,i),k=1,4)
299       enddo
300       nstrand=0
301       do i=1,nbfrag
302         write (iout,*) "calling add_strand:",i,bfrag(1,i),bfrag(2,i)
303         call add_strand(nstrand,istrand,nhairp,ihairp,
304      &     bfrag(1,i),bfrag(2,i),found)
305         if (bfrag(3,i).lt.bfrag(4,i)) then
306           write (iout,*) "calling add_strand:",i,bfrag(3,i),bfrag(4,i)
307           call add_strand(nstrand,istrand,nhairp,ihairp,
308      &     bfrag(3,i),bfrag(4,i),found)
309         else
310           write (iout,*) "calling add_strand:",i,bfrag(4,i),bfrag(3,i)
311           call add_strand(nstrand,istrand,nhairp,ihairp,
312      &      bfrag(4,i),bfrag(3,i),found)
313         endif
314       enddo
315       nbfrag=0
316       write (iout,*) "Strands found:",nstrand
317       do i=1,nstrand
318         write (iout,*) i,istrand(1,i),istrand(2,i)
319       enddo
320       return
321       end
322 c------------------------------------------------------------------------------
323       subroutine add_strand(nstrand,istrand,nhairp,ihairp,is1,is2,found)
324       implicit real*8 (a-h,o-z)
325       include 'DIMENSIONS'
326       include 'DIMENSIONS.ZSCOPT'
327       include 'DIMENSIONS.COMPAR'
328       include 'COMMON.IOUNITS'
329       integer nstrand,istrand(2,maxres/2)
330       integer nhairp,ihairp(2,maxres/5) 
331       logical found
332       found=.false.
333       do j=1,nhairp
334         idelt=(ihairp(2,j)-ihairp(1,j))/6
335         if (is1.lt.ihairp(2,j)-idelt.and.is2.gt.ihairp(1,j)+idelt) then
336           write (iout,*) "strand",is1,is2," is part of hairpin",
337      &      ihairp(1,j),ihairp(2,j)
338           return
339         endif
340       enddo
341       do j=1,nstrand
342         idelt=(istrand(2,j)-istrand(1,j))/3
343         if (is1.lt.istrand(2,j)-idelt.and.is2.gt.istrand(1,j)+idelt) 
344      &  then
345 c The strand already exists in the array; update its ends if necessary.
346           write (iout,*) "strand",is1,is2," found at position",j,
347      &     ":",istrand(1,j),istrand(2,j)
348           istrand(1,j)=min0(istrand(1,j),is1)
349           istrand(2,j)=max0(istrand(2,j),is2)
350           return   
351         endif
352       enddo
353 c The strand has not been found; add it to the array.
354       write (iout,*) "strand",is1,is2," added to the array."
355       found=.true.
356       nstrand=nstrand+1
357       istrand(1,nstrand)=is1
358       istrand(2,nstrand)=is2
359       return
360       end
361 c------------------------------------------------------------------------------
362       subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr)
363       implicit real*8 (a-h,o-z)
364       include 'DIMENSIONS'
365       include 'DIMENSIONS.ZSCOPT'
366       include 'COMMON.IOUNITS'
367       include 'COMMON.FRAG'
368       include 'COMMON.VAR'
369       include 'COMMON.GEO'
370       include 'COMMON.CHAIN'
371       include 'COMMON.NAMES'
372       include 'COMMON.INTERACT'
373       integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres),
374      &  isecstr(maxres)
375       logical lprint,lprint_sec,not_done,freeres
376       double precision p1,p2
377       external freeres
378       character*1 csec(0:2) /'-','E','H'/
379       if (lprint) then
380         write (iout,*) "entered secondary2",ncont
381         write (iout,*) "nstart_sup",nstart_sup," nend_sup",nend_sup
382         do i=1,ncont
383           write (iout,*) icont(1,i),icont(2,i)
384         enddo
385       endif
386       do i=1,nres
387         isecstr(i)=0
388       enddo
389       nbfrag=0
390       nhfrag=0
391       do i=1,nres
392         isec(i,1)=0
393         isec(i,2)=0
394         nsec(i)=0
395       enddo
396
397 c finding parallel beta
398 cd      write (iout,*) '------- looking for parallel beta -----------'
399       nbeta=0
400       nstrand=0
401       do i=1,ncont
402         i1=icont(1,i)
403         j1=icont(2,i)
404         if (i1.ge.nstart_sup .and. i1.le.nend_sup 
405      &     .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
406 cd        write (iout,*) "parallel",i1,j1
407         if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
408           ii1=i1
409           jj1=j1
410 cd          write (iout,*) i1,j1
411           not_done=.true.
412           do while (not_done)
413            i1=i1+1
414            j1=j1+1
415             do j=1,ncont
416               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
417      &             freeres(i1,j1,nsec,isec)) goto 5
418             enddo
419             not_done=.false.
420   5         continue
421 cd            write (iout,*) i1,j1,not_done
422           enddo
423           j1=j1-1
424           i1=i1-1
425           if (i1-ii1.gt.1) then
426             ii1=max0(ii1-1,1)
427             jj1=max0(jj1-1,1)
428             nbeta=nbeta+1
429             if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
430      &               nbeta,ii1,i1,jj1,j1
431
432             nbfrag=nbfrag+1
433             bfrag(1,nbfrag)=ii1+1
434             bfrag(2,nbfrag)=i1+1
435             bfrag(3,nbfrag)=jj1+1
436             bfrag(4,nbfrag)=min0(j1+1,nres) 
437
438             do ij=ii1,i1
439              nsec(ij)=nsec(ij)+1
440              isec(ij,nsec(ij))=nbeta
441             enddo
442             do ij=jj1,j1
443              nsec(ij)=nsec(ij)+1
444              isec(ij,nsec(ij))=nbeta
445             enddo
446
447            if(lprint_sec) then 
448             nstrand=nstrand+1
449             if (nbeta.le.9) then
450               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
451      &          "DefPropRes 'strand",nstrand,
452      &          "' 'num = ",ii1-1,"..",i1-1,"'"
453             else
454               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
455      &          "DefPropRes 'strand",nstrand,
456      &          "' 'num = ",ii1-1,"..",i1-1,"'"
457             endif
458             nstrand=nstrand+1
459             if (nbeta.le.9) then
460               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
461      &          "DefPropRes 'strand",nstrand,
462      &          "' 'num = ",jj1-1,"..",j1-1,"'"
463             else
464               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
465      &          "DefPropRes 'strand",nstrand,
466      &          "' 'num = ",jj1-1,"..",j1-1,"'"
467             endif
468               write(12,'(a8,4i4)')
469      &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
470            endif
471           endif
472         endif
473         endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup
474       enddo
475
476 c finding antiparallel beta
477 cd      write (iout,*) '--------- looking for antiparallel beta ---------'
478
479       do i=1,ncont
480         i1=icont(1,i)
481         j1=icont(2,i)
482         if (freeres(i1,j1,nsec,isec)) then
483           ii1=i1
484           jj1=j1
485 cd          write (iout,*) i1,j1
486
487           not_done=.true.
488           do while (not_done)
489            i1=i1+1
490            j1=j1-1
491             do j=1,ncont
492               if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
493      &             freeres(i1,j1,nsec,isec)) goto 6
494             enddo
495             not_done=.false.
496   6         continue
497 cd            write (iout,*) i1,j1,not_done
498           enddo
499           i1=i1-1
500           j1=j1+1
501           if (i1-ii1.gt.1) then
502
503             nbfrag=nbfrag+1
504             bfrag(1,nbfrag)=ii1
505             bfrag(2,nbfrag)=min0(i1+1,nres)
506             bfrag(3,nbfrag)=min0(jj1+1,nres)
507             bfrag(4,nbfrag)=j1
508
509             nbeta=nbeta+1
510             iii1=max0(ii1-1,1)
511             do ij=iii1,i1
512              nsec(ij)=nsec(ij)+1
513              if (nsec(ij).le.2) then
514               isec(ij,nsec(ij))=nbeta
515              endif
516             enddo
517             jjj1=max0(j1-1,1)  
518             do ij=jjj1,jj1
519              nsec(ij)=nsec(ij)+1
520              if (nsec(ij).le.2) then
521               isec(ij,nsec(ij))=nbeta
522              endif
523             enddo
524
525
526            if (lprint_sec) then
527             write (iout,'(a,i3,4i4)')'antiparallel beta',
528      &                   nbeta,ii1-1,i1,jj1,j1-1
529             nstrand=nstrand+1
530             if (nstrand.le.9) then
531               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
532      &          "DefPropRes 'strand",nstrand,
533      &          "' 'num = ",ii1-2,"..",i1-1,"'"
534             else
535               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
536      &          "DefPropRes 'strand",nstrand,
537      &          "' 'num = ",ii1-2,"..",i1-1,"'"
538             endif
539             nstrand=nstrand+1
540             if (nstrand.le.9) then
541               write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
542      &          "DefPropRes 'strand",nstrand,
543      &          "' 'num = ",j1-2,"..",jj1-1,"'"
544             else
545               write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
546      &          "DefPropRes 'strand",nstrand,
547      &          "' 'num = ",j1-2,"..",jj1-1,"'"
548             endif
549               write(12,'(a8,4i4)')
550      &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
551            endif
552           endif
553         endif
554       enddo
555
556 cd      write (iout,*) "After beta:",nbfrag
557 cd      do i=1,nbfrag
558 cd        write (iout,*) (bfrag(j,i),j=1,4)
559 cd      enddo
560
561       if (nstrand.gt.0.and.lprint_sec) then
562         write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
563         do i=2,nstrand
564          if (i.le.9) then
565           write(12,'(a9,i1,$)') " | strand",i
566          else
567           write(12,'(a9,i2,$)') " | strand",i
568          endif
569         enddo
570         write(12,'(a1)') "'"
571       endif
572
573        
574 c finding alpha or 310 helix
575
576       nhelix=0
577       do i=1,ncont
578         i1=icont(1,i)
579         j1=icont(2,i)
580         p1=phi(i1+2)*rad2deg
581         p2=0.0
582         if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
583
584
585         if (j1.eq.i1+3 .and. 
586      &       ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
587      &       ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
588 cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
589 co          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
590           ii1=i1
591           jj1=j1
592           if (nsec(ii1).eq.0) then 
593             not_done=.true.
594           else
595             not_done=.false.
596           endif
597           do while (not_done)
598             i1=i1+1
599             j1=j1+1
600             do j=1,ncont
601               if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
602             enddo
603             not_done=.false.
604   10        continue
605             p1=phi(i1+2)*rad2deg
606             p2=phi(j1+2)*rad2deg
607             if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) 
608      &                              not_done=.false.
609
610 cd           write (iout,*) i1,j1,not_done,p1,p2
611           enddo
612           j1=j1+1
613           if (j1-ii1.gt.4) then
614             nhelix=nhelix+1
615 cd            write (iout,*)'helix',nhelix,ii1,j1
616
617             nhfrag=nhfrag+1
618             hfrag(1,nhfrag)=ii1
619             hfrag(2,nhfrag)=j1
620
621             do ij=ii1,j1
622              nsec(ij)=-1
623             enddo
624            if (lprint_sec) then
625             write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
626             if (nhelix.le.9) then
627               write(12,'(a17,i1,a9,i3,a2,i3,a1)') 
628      &          "DefPropRes 'helix",nhelix,
629      &          "' 'num = ",ii1-1,"..",j1-2,"'"
630             else
631               write(12,'(a17,i2,a9,i3,a2,i3,a1)') 
632      &          "DefPropRes 'helix",nhelix,
633      &          "' 'num = ",ii1-1,"..",j1-2,"'"
634             endif
635            endif
636           endif
637         endif
638       enddo
639        
640       if (nhelix.gt.0.and.lprint_sec) then
641         write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
642         do i=2,nhelix
643          if (nhelix.le.9) then
644           write(12,'(a8,i1,$)') " | helix",i
645          else
646           write(12,'(a8,i2,$)') " | helix",i
647          endif
648         enddo
649         write(12,'(a1)') "'"
650       endif
651
652       if (lprint_sec) then
653        write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
654        write(12,'(a20)') "XMacStand ribbon.mac"
655       endif
656         
657       if (lprint) then
658
659         write(iout,*) 'UNRES seq:'
660         do j=1,nbfrag
661          write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
662         enddo
663   
664         do j=1,nhfrag
665          write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
666         enddo
667
668       endif   
669   
670       do j=1,nbfrag
671         do k=min0(bfrag(1,j),bfrag(2,j)),max0(bfrag(1,j),bfrag(2,j)) 
672           isecstr(k)=1
673         enddo
674         do k=min0(bfrag(3,j),bfrag(4,j)),max0(bfrag(3,j),bfrag(4,j)) 
675           isecstr(k)=1
676         enddo
677       enddo
678       do j=1,nhfrag
679         do k=hfrag(1,j),hfrag(2,j)
680           isecstr(k)=2
681         enddo
682       enddo
683       if (lprint) then
684         write (iout,*)
685         write (iout,*) "Secondary structure"
686         do i=1,nres,80
687           ist=i
688           ien=min0(i+79,nres)
689           write (iout,*)
690           write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10)
691           write (iout,'(80a1)') (onelet(itype(k)),k=ist,ien) 
692           write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien) 
693         enddo 
694         write (iout,*)
695       endif
696       return
697       end
698 c-------------------------------------------------
699       logical function freeres(i,j,nsec,isec)
700       include 'DIMENSIONS'
701       integer isec(maxres,4),nsec(maxres)
702       freeres=.false.
703
704       if (nsec(i).gt.1.or.nsec(j).gt.1) return
705       do k=1,nsec(i)
706         do l=1,nsec(j)
707           if (isec(i,k).eq.isec(j,l)) return
708         enddo
709       enddo
710       freeres=.true.
711       return
712       end
713