X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-NEWSC%2Fsecondary.f;fp=source%2Fwham%2Fsrc-NEWSC%2Fsecondary.f;h=9c9bc7d949d9c191d1c7f62ca31b143fa28e8c56;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git

diff --git a/source/wham/src-NEWSC/secondary.f b/source/wham/src-NEWSC/secondary.f
new file mode 100755
index 0000000..9c9bc7d
--- /dev/null
+++ b/source/wham/src-NEWSC/secondary.f
@@ -0,0 +1,713 @@
+      subroutine define_fragments
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.COMPAR'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.TIME1'
+      include 'COMMON.FRAG'
+      include 'COMMON.SBRIDGE'
+      include 'COMMON.CONTROL'
+      include 'COMMON.COMPAR'
+      include 'COMMON.CHAIN'
+      include 'COMMON.HEADER'
+      include 'COMMON.GEO'
+      include 'COMMON.CONTACTS'
+      include 'COMMON.PEPTCONT'
+      include 'COMMON.INTERACT'
+      include 'COMMON.NAMES'
+      integer nstrand,istrand(2,maxres/2)
+      integer nhairp,ihairp(2,maxres/5) 
+      character*16 strstr(4) /'helix','hairpin','strand','strand pair'/
+      write (iout,*) 'NC_FRAC_HEL',ncfrac_hel,' NC_REQ_HEL',ncreq_hel,
+     &               'NC_FRAC_BET',ncfrac_bet,' NC_REQ_BET',ncreq_bet,
+     &           'NC_FRAC_PAIR',ncfrac_pair,' NC_REQ_PAIR',ncreq_pair,
+     &  ' RMS_PAIR',irms_pair,' SPLIT_BET',isplit_bet
+      write (iout,*) 'NSHIFT_HEL',nshift_hel,' NSHIFT_BET',nshift_bet,
+     &  ' NSHIFT_STRAND',nshift_strand,' NSHIFT_PAIR',nshift_pair
+      write (iout,*) 'ANGCUT_HEL',angcut_hel*rad2deg,
+     &  ' MAXANG_HEL',angcut1_hel*rad2deg
+      write (iout,*) 'ANGCUT_BET',angcut_bet*rad2deg,
+     &               ' MAXANG_BET',angcut1_bet*rad2deg
+      write (iout,*) 'ANGCUT_STRAND',angcut_strand*rad2deg,
+     &               ' MAXANG_STRAND',angcut1_strand*rad2deg
+      write (iout,*) 'FRAC_MIN',frac_min_set
+c Find secondary structure elements (helices and beta-sheets)
+      call secondary2(.true.,.false.,ncont_pept_ref,icont_pept_ref,
+     &   isec_ref)
+c Define primary fragments. First include the helices.
+      nhairp=0
+      nstrand=0
+c Merge helices
+c AL 12/23/03 - to avoid splitting helices into very small fragments
+      if (merge_helices) then
+      write (iout,*) "Before merging helices: nhfrag",nhfrag
+      do i=1,nhfrag
+        write (2,*) hfrag(1,i),hfrag(2,i)
+      enddo
+      i=1
+      do while (i.lt.nhfrag)
+        if (hfrag(1,i+1)-hfrag(2,i).le.1) then
+          nhfrag=nhfrag-1
+          hfrag(2,i)=hfrag(2,i+1)
+          do j=i+1,nhfrag
+            hfrag(1,j)=hfrag(1,j+1)
+            hfrag(2,j)=hfrag(2,j+1)
+          enddo
+        endif 
+        i=i+1
+      enddo
+      write (iout,*) "After merging helices: nhfrag",nhfrag
+      do i=1,nhfrag
+        write (2,*) hfrag(1,i),hfrag(2,i)
+      enddo
+      endif
+      nfrag(1)=nhfrag
+      do i=1,nhfrag
+        npiece(i,1)=1
+        ifrag(1,1,i)=hfrag(1,i) 
+        ifrag(2,1,i)=hfrag(2,i) 
+        n_shift(1,i,1)=0
+        n_shift(2,i,1)=nshift_hel
+        ang_cut(i)=angcut_hel
+        ang_cut1(i)=angcut1_hel
+        frac_min(i)=frac_min_set
+        nc_fragm(i,1)=ncfrac_hel
+        nc_req_setf(i,1)=ncreq_hel
+        istruct(i)=1
+      enddo
+      write (iout,*) "isplit_bet",isplit_bet
+      if (isplit_bet.gt.1) then
+c Split beta-sheets into strands and store strands as primary fragments.
+        call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
+        do i=1,nstrand
+          ii=i+nfrag(1)
+          npiece(ii,1)=1
+          ifrag(1,1,ii)=istrand(1,i)
+          ifrag(2,1,ii)=istrand(2,i)
+          n_shift(1,ii,1)=nshift_strand
+          n_shift(2,ii,1)=nshift_strand
+          ang_cut(ii)=angcut_strand
+          ang_cut1(ii)=angcut1_strand
+          frac_min(ii)=frac_min_set
+          nc_fragm(ii,1)=0
+          nc_req_setf(ii,1)=0
+          istruct(ii)=3
+        enddo
+        nfrag(1)=nfrag(1)+nstrand
+      else if (isplit_bet.eq.1) then
+c Split only far beta-sheets; does not split hairpins.
+        call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
+        call split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
+        do i=1,nhairp
+          ii=i+nfrag(1)
+          npiece(ii,1)=1
+          ifrag(1,1,ii)=ihairp(1,i) 
+          ifrag(2,1,ii)=ihairp(2,i) 
+          n_shift(1,ii,1)=nshift_bet
+          n_shift(2,ii,1)=nshift_bet
+          ang_cut(ii)=angcut_bet
+          ang_cut1(ii)=angcut1_bet
+          frac_min(ii)=frac_min_set
+          nc_fragm(ii,1)=ncfrac_bet
+          nc_req_setf(ii,1)=ncreq_bet
+          istruct(ii)=2
+        enddo
+        nfrag(1)=nfrag(1)+nhairp
+        do i=1,nstrand
+          ii=i+nfrag(1)
+          npiece(ii,1)=1
+          ifrag(1,1,ii)=istrand(1,i)
+          ifrag(2,1,ii)=istrand(2,i)
+          n_shift(1,ii,1)=nshift_strand
+          n_shift(2,ii,1)=nshift_strand
+          ang_cut(ii)=angcut_strand
+          ang_cut1(ii)=angcut1_strand
+          frac_min(ii)=frac_min_set
+          nc_fragm(ii,1)=0
+          nc_req_setf(ii,1)=0
+          istruct(ii)=3
+        enddo
+        nfrag(1)=nfrag(1)+nstrand
+      else
+c Do not split beta-sheets; each pair of strands is a primary element.
+        call find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
+        do i=1,nhairp
+          ii=i+nfrag(1)
+          npiece(ii,1)=1
+          ifrag(1,1,ii)=ihairp(1,i) 
+          ifrag(2,1,ii)=ihairp(2,i) 
+          n_shift(1,ii,1)=nshift_bet
+          n_shift(2,ii,1)=nshift_bet
+          ang_cut(ii)=angcut_bet
+          ang_cut1(ii)=angcut1_bet
+          frac_min(ii)=frac_min_set
+          nc_fragm(ii,1)=ncfrac_bet
+          nc_req_setf(ii,1)=ncreq_bet
+          istruct(ii)=2
+        enddo
+        nfrag(1)=nfrag(1)+nhairp
+        do i=1,nbfrag
+          ii=i+nfrag(1)
+          npiece(ii,1)=2
+          ifrag(1,1,ii)=bfrag(1,i) 
+          ifrag(2,1,ii)=bfrag(2,i) 
+          if (bfrag(3,i).lt.bfrag(4,i)) then
+            ifrag(1,2,ii)=bfrag(3,i)
+            ifrag(2,2,ii)=bfrag(4,i)
+          else
+            ifrag(1,2,ii)=bfrag(4,i)
+            ifrag(2,2,ii)=bfrag(3,i)
+          endif
+          n_shift(1,ii,1)=nshift_bet
+          n_shift(2,ii,1)=nshift_bet
+          ang_cut(ii)=angcut_bet
+          ang_cut1(ii)=angcut1_bet
+          frac_min(ii)=frac_min_set
+          nc_fragm(ii,1)=ncfrac_bet
+          nc_req_setf(ii,1)=ncreq_bet
+          istruct(ii)=4
+        enddo
+        nfrag(1)=nfrag(1)+nbfrag
+      endif
+      write (iout,*) "The following primary fragments were found:"
+      write (iout,*) "Helices:",nhfrag
+      do i=1,nhfrag
+        i1=ifrag(1,1,i)
+        i2=ifrag(2,1,i)
+        it1=itype(i1)
+        it2=itype(i2)
+        write (iout,'(i3,2x,a,i4,2x,a,i4)')
+     &       i,restyp(it1),i1,restyp(it2),i2
+      enddo
+      write (iout,*) "Hairpins:",nhairp
+      do i=nhfrag+1,nhfrag+nhairp
+        i1=ifrag(1,1,i)
+        i2=ifrag(2,1,i)
+        it1=itype(i1)
+        it2=itype(i2)
+        write (iout,'(i3,2x,a,i4,2x,a,i4,2x)')
+     &       i,restyp(it1),i1,restyp(it2),i2
+      enddo
+      write (iout,*) "Far strand pairs:",nbfrag
+      do i=nhfrag+nhairp+1,nhfrag+nhairp+nbfrag
+        i1=ifrag(1,1,i)
+        i2=ifrag(2,1,i)
+        it1=itype(i1)
+        it2=itype(i2)
+        i3=ifrag(1,2,i)
+        i4=ifrag(2,2,i)
+        it3=itype(i3)
+        it4=itype(i4)
+        write (iout,'(i3,2x,a,i4,2x,a,i4," and ",a,i4,2x,a,i4)')
+     &       i,restyp(it1),i1,restyp(it2),i2,
+     &         restyp(it3),i3,restyp(it4),i4
+      enddo
+      write (iout,*) "Strands:",nstrand
+      do i=nhfrag+nhairp+nbfrag+1,nfrag(1)
+        i1=ifrag(1,1,i)
+        i2=ifrag(2,1,i)
+        it1=itype(i1)
+        it2=itype(i2)
+        write (iout,'(i3,2x,a,i4,2x,a,i4)')
+     &       i,restyp(it1),i1,restyp(it2),i2
+      enddo
+      call imysort(nfrag(1),2,maxpiece,ifrag(1,1,1),npiece(1,1),
+     &  istruct(1),n_shift(1,1,1),ang_cut(1),ang_cut1(1),frac_min(1),
+     &  nc_fragm(1,1),nc_req_setf(1,1))
+      write (iout,*) "Fragments after sorting:"
+      do i=1,nfrag(1)
+        i1=ifrag(1,1,i)
+        i2=ifrag(2,1,i)
+        it1=itype(i1)
+        it2=itype(i2)
+        write (iout,'(i3,2x,a,i4,2x,a,i4,$)')
+     &       i,restyp(it1),i1,restyp(it2),i2
+        if (npiece(i,1).eq.1) then
+          write (iout,'(2x,a)') strstr(istruct(i))
+        else
+          i1=ifrag(1,2,i)
+          i2=ifrag(2,2,i)
+          it1=itype(i1)
+          it2=itype(i2)
+          write (iout,'(2x,a,i4,2x,a,i4,2x,a)')
+     &       restyp(it1),i1,restyp(it2),i2,strstr(istruct(i))
+        endif
+      enddo
+      return
+      end
+c------------------------------------------------------------------------------
+      subroutine find_and_remove_hairpins(nbfrag,bfrag,nhairp,ihairp)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.COMPAR'
+      include 'COMMON.IOUNITS'
+      integer nbfrag,bfrag(4,maxres/3)
+      integer nhairp,ihairp(2,maxres/5) 
+      write (iout,*) "Entered find_and_remove_hairpins"
+      write (iout,*) "nbfrag",nbfrag
+      do i=1,nbfrag
+        write (iout,*) i,(bfrag(k,i),k=1,4)
+      enddo
+      nhairp=0
+      i=1
+      do while (i.le.nbfrag)
+        write (iout,*) "check hairpin:",i,(bfrag(j,i),j=1,4)
+        if (bfrag(3,i).gt.bfrag(4,i) .and. bfrag(4,i)-bfrag(2,i).lt.5) 
+     &  then
+          write (iout,*) "Found hairpin:",i,bfrag(1,i),bfrag(3,i)
+          nhairp=nhairp+1
+          ihairp(1,nhairp)=bfrag(1,i)
+          ihairp(2,nhairp)=bfrag(3,i) 
+          nbfrag=nbfrag-1
+          do j=i,nbfrag
+            do k=1,4
+              bfrag(k,j)=bfrag(k,j+1)
+            enddo
+          enddo
+        else
+          i=i+1
+        endif
+      enddo
+      write (iout,*) "After finding hairpins:"
+      write (iout,*) "nhairp",nhairp
+      do i=1,nhairp
+        write (iout,*) i,ihairp(1,i),ihairp(2,i)
+      enddo
+      write (iout,*) "nbfrag",nbfrag
+      do i=1,nbfrag
+        write (iout,*) i,(bfrag(k,i),k=1,4)
+      enddo
+      return
+      end
+c------------------------------------------------------------------------------
+      subroutine split_beta(nbfrag,bfrag,nstrand,istrand,nhairp,ihairp)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.COMPAR'
+      include 'COMMON.IOUNITS'
+      integer nbfrag,bfrag(4,maxres/3)
+      integer nstrand,istrand(2,maxres/2)
+      integer nhairp,ihairp(2,maxres/5) 
+      logical found
+      write (iout,*) "Entered split_beta"
+      write (iout,*) "nbfrag",nbfrag
+      do i=1,nbfrag
+        write (iout,*) i,(bfrag(k,i),k=1,4)
+      enddo
+      nstrand=0
+      do i=1,nbfrag
+        write (iout,*) "calling add_strand:",i,bfrag(1,i),bfrag(2,i)
+        call add_strand(nstrand,istrand,nhairp,ihairp,
+     &     bfrag(1,i),bfrag(2,i),found)
+        if (bfrag(3,i).lt.bfrag(4,i)) then
+          write (iout,*) "calling add_strand:",i,bfrag(3,i),bfrag(4,i)
+          call add_strand(nstrand,istrand,nhairp,ihairp,
+     &     bfrag(3,i),bfrag(4,i),found)
+        else
+          write (iout,*) "calling add_strand:",i,bfrag(4,i),bfrag(3,i)
+          call add_strand(nstrand,istrand,nhairp,ihairp,
+     &      bfrag(4,i),bfrag(3,i),found)
+        endif
+      enddo
+      nbfrag=0
+      write (iout,*) "Strands found:",nstrand
+      do i=1,nstrand
+        write (iout,*) i,istrand(1,i),istrand(2,i)
+      enddo
+      return
+      end
+c------------------------------------------------------------------------------
+      subroutine add_strand(nstrand,istrand,nhairp,ihairp,is1,is2,found)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'DIMENSIONS.COMPAR'
+      include 'COMMON.IOUNITS'
+      integer nstrand,istrand(2,maxres/2)
+      integer nhairp,ihairp(2,maxres/5) 
+      logical found
+      found=.false.
+      do j=1,nhairp
+        idelt=(ihairp(2,j)-ihairp(1,j))/6
+        if (is1.lt.ihairp(2,j)-idelt.and.is2.gt.ihairp(1,j)+idelt) then
+          write (iout,*) "strand",is1,is2," is part of hairpin",
+     &      ihairp(1,j),ihairp(2,j)
+          return
+        endif
+      enddo
+      do j=1,nstrand
+        idelt=(istrand(2,j)-istrand(1,j))/3
+        if (is1.lt.istrand(2,j)-idelt.and.is2.gt.istrand(1,j)+idelt) 
+     &  then
+c The strand already exists in the array; update its ends if necessary.
+          write (iout,*) "strand",is1,is2," found at position",j,
+     &     ":",istrand(1,j),istrand(2,j)
+          istrand(1,j)=min0(istrand(1,j),is1)
+          istrand(2,j)=max0(istrand(2,j),is2)
+          return   
+        endif
+      enddo
+c The strand has not been found; add it to the array.
+      write (iout,*) "strand",is1,is2," added to the array."
+      found=.true.
+      nstrand=nstrand+1
+      istrand(1,nstrand)=is1
+      istrand(2,nstrand)=is2
+      return
+      end
+c------------------------------------------------------------------------------
+      subroutine secondary2(lprint,lprint_sec,ncont,icont,isecstr)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'DIMENSIONS.ZSCOPT'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.FRAG'
+      include 'COMMON.VAR'
+      include 'COMMON.GEO'
+      include 'COMMON.CHAIN'
+      include 'COMMON.NAMES'
+      include 'COMMON.INTERACT'
+      integer ncont,icont(2,maxcont),isec(maxres,4),nsec(maxres),
+     &  isecstr(maxres)
+      logical lprint,lprint_sec,not_done,freeres
+      double precision p1,p2
+      external freeres
+      character*1 csec(0:2) /'-','E','H'/
+      if (lprint) then
+        write (iout,*) "entered secondary2",ncont
+        write (iout,*) "nstart_sup",nstart_sup," nend_sup",nend_sup
+        do i=1,ncont
+          write (iout,*) icont(1,i),icont(2,i)
+        enddo
+      endif
+      do i=1,nres
+        isecstr(i)=0
+      enddo
+      nbfrag=0
+      nhfrag=0
+      do i=1,nres
+        isec(i,1)=0
+        isec(i,2)=0
+        nsec(i)=0
+      enddo
+
+c finding parallel beta
+cd      write (iout,*) '------- looking for parallel beta -----------'
+      nbeta=0
+      nstrand=0
+      do i=1,ncont
+        i1=icont(1,i)
+        j1=icont(2,i)
+        if (i1.ge.nstart_sup .and. i1.le.nend_sup 
+     &     .and. j1.gt.nstart_sup .and. j1.le.nend_sup) then
+cd        write (iout,*) "parallel",i1,j1
+        if(j1-i1.gt.5 .and. freeres(i1,j1,nsec,isec)) then
+          ii1=i1
+          jj1=j1
+cd          write (iout,*) i1,j1
+          not_done=.true.
+          do while (not_done)
+           i1=i1+1
+           j1=j1+1
+            do j=1,ncont
+              if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j) .and.
+     &             freeres(i1,j1,nsec,isec)) goto 5
+            enddo
+            not_done=.false.
+  5         continue
+cd            write (iout,*) i1,j1,not_done
+          enddo
+          j1=j1-1
+          i1=i1-1
+          if (i1-ii1.gt.1) then
+            ii1=max0(ii1-1,1)
+            jj1=max0(jj1-1,1)
+            nbeta=nbeta+1
+            if(lprint)write(iout,'(a,i3,4i4)')'parallel beta',
+     &               nbeta,ii1,i1,jj1,j1
+
+            nbfrag=nbfrag+1
+            bfrag(1,nbfrag)=ii1+1
+            bfrag(2,nbfrag)=i1+1
+            bfrag(3,nbfrag)=jj1+1
+            bfrag(4,nbfrag)=min0(j1+1,nres) 
+
+            do ij=ii1,i1
+             nsec(ij)=nsec(ij)+1
+             isec(ij,nsec(ij))=nbeta
+            enddo
+            do ij=jj1,j1
+             nsec(ij)=nsec(ij)+1
+             isec(ij,nsec(ij))=nbeta
+            enddo
+
+           if(lprint_sec) then 
+            nstrand=nstrand+1
+            if (nbeta.le.9) then
+              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",ii1-1,"..",i1-1,"'"
+            else
+              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",ii1-1,"..",i1-1,"'"
+            endif
+            nstrand=nstrand+1
+            if (nbeta.le.9) then
+              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",jj1-1,"..",j1-1,"'"
+            else
+              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",jj1-1,"..",j1-1,"'"
+            endif
+              write(12,'(a8,4i4)')
+     &          "SetNeigh",ii1-1,i1-1,jj1-1,j1-1
+           endif
+          endif
+        endif
+        endif ! i1.ge.nstart_sup .and. i1.le.nend_sup .and. i2.gt.nstart_sup .and. i2.le.nend_sup
+      enddo
+
+c finding antiparallel beta
+cd      write (iout,*) '--------- looking for antiparallel beta ---------'
+
+      do i=1,ncont
+        i1=icont(1,i)
+        j1=icont(2,i)
+        if (freeres(i1,j1,nsec,isec)) then
+          ii1=i1
+          jj1=j1
+cd          write (iout,*) i1,j1
+
+          not_done=.true.
+          do while (not_done)
+           i1=i1+1
+           j1=j1-1
+            do j=1,ncont
+              if (i1.eq.icont(1,j).and.j1.eq.icont(2,j) .and.
+     &             freeres(i1,j1,nsec,isec)) goto 6
+            enddo
+            not_done=.false.
+  6         continue
+cd            write (iout,*) i1,j1,not_done
+          enddo
+          i1=i1-1
+          j1=j1+1
+          if (i1-ii1.gt.1) then
+
+            nbfrag=nbfrag+1
+            bfrag(1,nbfrag)=ii1
+            bfrag(2,nbfrag)=min0(i1+1,nres)
+            bfrag(3,nbfrag)=min0(jj1+1,nres)
+            bfrag(4,nbfrag)=j1
+
+            nbeta=nbeta+1
+            iii1=max0(ii1-1,1)
+            do ij=iii1,i1
+             nsec(ij)=nsec(ij)+1
+             if (nsec(ij).le.2) then
+              isec(ij,nsec(ij))=nbeta
+             endif
+            enddo
+            jjj1=max0(j1-1,1)  
+            do ij=jjj1,jj1
+             nsec(ij)=nsec(ij)+1
+             if (nsec(ij).le.2) then
+              isec(ij,nsec(ij))=nbeta
+             endif
+            enddo
+
+
+           if (lprint_sec) then
+            write (iout,'(a,i3,4i4)')'antiparallel beta',
+     &                   nbeta,ii1-1,i1,jj1,j1-1
+            nstrand=nstrand+1
+            if (nstrand.le.9) then
+              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",ii1-2,"..",i1-1,"'"
+            else
+              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",ii1-2,"..",i1-1,"'"
+            endif
+            nstrand=nstrand+1
+            if (nstrand.le.9) then
+              write(12,'(a18,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",j1-2,"..",jj1-1,"'"
+            else
+              write(12,'(a18,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'strand",nstrand,
+     &          "' 'num = ",j1-2,"..",jj1-1,"'"
+            endif
+              write(12,'(a8,4i4)')
+     &          "SetNeigh",ii1-2,i1-1,jj1-1,j1-2
+           endif
+          endif
+        endif
+      enddo
+
+cd      write (iout,*) "After beta:",nbfrag
+cd      do i=1,nbfrag
+cd        write (iout,*) (bfrag(j,i),j=1,4)
+cd      enddo
+
+      if (nstrand.gt.0.and.lprint_sec) then
+        write(12,'(a27,$)') "DefPropRes 'sheet' 'strand1"
+        do i=2,nstrand
+         if (i.le.9) then
+          write(12,'(a9,i1,$)') " | strand",i
+         else
+          write(12,'(a9,i2,$)') " | strand",i
+         endif
+        enddo
+        write(12,'(a1)') "'"
+      endif
+
+       
+c finding alpha or 310 helix
+
+      nhelix=0
+      do i=1,ncont
+        i1=icont(1,i)
+        j1=icont(2,i)
+        p1=phi(i1+2)*rad2deg
+        p2=0.0
+        if (j1+2.le.nres) p2=phi(j1+2)*rad2deg
+
+
+        if (j1.eq.i1+3 .and. 
+     &       ((p1.ge.10.and.p1.le.80).or.i1.le.2).and.
+     &       ((p2.ge.10.and.p2.le.80).or.j1.le.2.or.j1.ge.nres-3) )then
+cd          if (j1.eq.i1+3) write (iout,*) "found 1-4 ",i1,j1,p1,p2
+co          if (j1.eq.i1+4) write (iout,*) "found 1-5 ",i1,j1,p1,p2
+          ii1=i1
+          jj1=j1
+          if (nsec(ii1).eq.0) then 
+            not_done=.true.
+          else
+            not_done=.false.
+          endif
+          do while (not_done)
+            i1=i1+1
+            j1=j1+1
+            do j=1,ncont
+              if (i1.eq.icont(1,j) .and. j1.eq.icont(2,j)) goto 10
+            enddo
+            not_done=.false.
+  10        continue
+            p1=phi(i1+2)*rad2deg
+            p2=phi(j1+2)*rad2deg
+            if (p1.lt.10.or.p1.gt.80.or.p2.lt.10.or.p2.gt.80) 
+     &                              not_done=.false.
+
+cd           write (iout,*) i1,j1,not_done,p1,p2
+          enddo
+          j1=j1+1
+          if (j1-ii1.gt.4) then
+            nhelix=nhelix+1
+cd            write (iout,*)'helix',nhelix,ii1,j1
+
+            nhfrag=nhfrag+1
+            hfrag(1,nhfrag)=ii1
+            hfrag(2,nhfrag)=j1
+
+            do ij=ii1,j1
+             nsec(ij)=-1
+            enddo
+           if (lprint_sec) then
+            write (iout,'(a,i3,2i4)') "Helix",nhelix,ii1-1,j1-1
+            if (nhelix.le.9) then
+              write(12,'(a17,i1,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'helix",nhelix,
+     &          "' 'num = ",ii1-1,"..",j1-2,"'"
+            else
+              write(12,'(a17,i2,a9,i3,a2,i3,a1)') 
+     &          "DefPropRes 'helix",nhelix,
+     &          "' 'num = ",ii1-1,"..",j1-2,"'"
+            endif
+           endif
+          endif
+        endif
+      enddo
+       
+      if (nhelix.gt.0.and.lprint_sec) then
+        write(12,'(a26,$)') "DefPropRes 'helix' 'helix1"
+        do i=2,nhelix
+         if (nhelix.le.9) then
+          write(12,'(a8,i1,$)') " | helix",i
+         else
+          write(12,'(a8,i2,$)') " | helix",i
+         endif
+        enddo
+        write(12,'(a1)') "'"
+      endif
+
+      if (lprint_sec) then
+       write(12,'(a37)') "DefPropRes 'coil' '! (helix | sheet)'"
+       write(12,'(a20)') "XMacStand ribbon.mac"
+      endif
+        
+      if (lprint) then
+
+        write(iout,*) 'UNRES seq:'
+        do j=1,nbfrag
+         write(iout,*) 'beta ',(bfrag(i,j),i=1,4)
+        enddo
+  
+        do j=1,nhfrag
+         write(iout,*) 'helix ',(hfrag(i,j),i=1,2)
+        enddo
+
+      endif   
+  
+      do j=1,nbfrag
+        do k=min0(bfrag(1,j),bfrag(2,j)),max0(bfrag(1,j),bfrag(2,j)) 
+          isecstr(k)=1
+        enddo
+        do k=min0(bfrag(3,j),bfrag(4,j)),max0(bfrag(3,j),bfrag(4,j)) 
+          isecstr(k)=1
+        enddo
+      enddo
+      do j=1,nhfrag
+        do k=hfrag(1,j),hfrag(2,j)
+          isecstr(k)=2
+        enddo
+      enddo
+      if (lprint) then
+        write (iout,*)
+        write (iout,*) "Secondary structure"
+        do i=1,nres,80
+          ist=i
+          ien=min0(i+79,nres)
+          write (iout,*)
+          write (iout,'(8(7x,i3))') (k,k=ist+9,ien,10)
+          write (iout,'(80a1)') (onelet(itype(k)),k=ist,ien) 
+          write (iout,'(80a1)') (csec(isecstr(k)),k=ist,ien) 
+        enddo 
+        write (iout,*)
+      endif
+      return
+      end
+c-------------------------------------------------
+      logical function freeres(i,j,nsec,isec)
+      include 'DIMENSIONS'
+      integer isec(maxres,4),nsec(maxres)
+      freeres=.false.
+
+      if (nsec(i).gt.1.or.nsec(j).gt.1) return
+      do k=1,nsec(i)
+        do l=1,nsec(j)
+          if (isec(i,k).eq.isec(j,l)) return
+        enddo
+      enddo
+      freeres=.true.
+      return
+      end
+