source/wham/src-M-Pawel/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:',anatemp
 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),anatemp
 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