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