source/wham/src-M-homology/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 'DIMENSIONS.FREE'
   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