+
+!-----------------------------------------------------------
+ subroutine seq2chains
+!c
+!c Split the total UNRES sequence, which has dummy residues separating
+!c the chains, into separate chains. The length of chain ichain is
+!c contained in chain_length(ichain), the first and last non-dummy
+!c residues are in chain_border(1,ichain) and chain_border(2,ichain),
+!c respectively. The lengths pertain to non-dummy residues only.
+!c
+! implicit none
+! include 'DIMENSIONS'
+ use energy_data, only:molnum,nchain,chain_length,ireschain
+ implicit none
+! integer ireschain(nres)
+ integer ii,ichain,i,j,mnum
+ logical new_chain
+ print *,"in seq2"
+ ichain=1
+ new_chain=.true.
+ if (.not.allocated(chain_length)) allocate(chain_length(50))
+ if (.not.allocated(chain_border)) allocate(chain_border(2,50))
+
+ chain_length(ichain)=0
+ ii=1
+ do while (ii.lt.nres)
+ write(iout,*) "in seq2chains",ii,new_chain
+ mnum=molnum(ii)
+ if (itype(ii,mnum).eq.ntyp1_molec(mnum)) then
+ if (.not.new_chain) then
+ new_chain=.true.
+ chain_border(2,ichain)=ii-1
+ ichain=ichain+1
+ chain_border(1,ichain)=ii+1
+ chain_length(ichain)=0
+ endif
+ else
+ if (new_chain) then
+ chain_border(1,ichain)=ii
+ new_chain=.false.
+ endif
+ chain_length(ichain)=chain_length(ichain)+1
+ endif
+ ii=ii+1
+ enddo
+ if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) then
+ ii=ii-1
+ else
+ chain_length(ichain)=chain_length(ichain)+1
+ endif
+ if (chain_length(ichain).gt.0) then
+ chain_border(2,ichain)=ii
+ nchain=ichain
+ else
+ nchain=ichain-1
+ endif
+ ireschain=0
+ do i=1,nchain
+ do j=chain_border(1,i),chain_border(2,i)
+ ireschain(j)=i
+ enddo
+ enddo
+ return
+ end subroutine
+!---------------------------------------------------------------------
+ subroutine chain_symmetry(npermchain,tabpermchain)
+!c
+!c Determine chain symmetry. nperm is the number of permutations and
+!c tabperchain contains the allowed permutations of the chains.
+!c
+! implicit none
+! include "DIMENSIONS"
+! include "COMMON.IOUNITS"
+ implicit none
+ integer itemp(50),&
+ npermchain,tabpermchain(50,5040),&
+ tabperm(50,5040),mapchain(50),&
+ iequiv(50,nres),iflag(nres)
+ integer i,j,k,l,ii,nchain_group,nequiv(50),iieq,&
+ nperm,npermc,ind,mnum
+ if (nchain.eq.1) then
+ npermchain=1
+ tabpermchain(1,1)=1
+!c print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
+ return
+ endif
+!c
+!c Look for equivalent chains
+#ifdef DEBUG
+ write(iout,*) "nchain",nchain
+ do i=1,nchain
+ write(iout,*) "chain",i," from",chain_border(1,i),&
+ " to",chain_border(2,i)
+ write(iout,*)&
+ "sequence ",(itype(j,molnum(j)),j=chain_border(1,i),chain_border(2,i))
+ enddo
+#endif
+ do i=1,nchain
+ iflag(i)=0
+ enddo
+ nchain_group=0
+ do i=1,nchain
+ if (iflag(i).gt.0) cycle
+ iflag(i)=1
+ nchain_group=nchain_group+1
+ iieq=1
+ iequiv(iieq,nchain_group)=i
+ do j=i+1,nchain
+ if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
+!c k=0
+!c do while(k.lt.chain_length(i) .and.
+!c & itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
+ do k=0,chain_length(i)-1
+!c k=k+1
+ mnum=molnum(k)
+ if (itype(chain_border(1,i)+k,mnum).ne.&
+ itype(chain_border(1,j)+k,mnum)) exit
+ enddo
+ if (k.lt.chain_length(i)) cycle
+ iflag(j)=1
+ iieq=iieq+1
+ iequiv(iieq,nchain_group)=j
+ enddo
+ nequiv(nchain_group)=iieq
+ enddo
+ write(iout,*) "Number of equivalent chain groups:",nchain_group
+ write(iout,*) "Equivalent chain groups"
+ do i=1,nchain_group
+ write(iout,*) "group",i," #members",nequiv(i)," chains",&
+ (iequiv(j,i),j=1,nequiv(i))
+ enddo
+ ind=0
+ do i=1,nchain_group
+ do j=1,nequiv(i)
+ ind=ind+1
+ mapchain(ind)=iequiv(j,i)
+ enddo
+ enddo
+ write (iout,*) "mapchain"
+ do i=1,nchain
+ write (iout,*) i,mapchain(i)
+ enddo
+ ii=0
+ do i=1,nchain_group
+ call permut(nequiv(i),nperm,tabperm)
+ if (ii.eq.0) then
+ ii=nequiv(i)
+ npermchain=nperm
+ do j=1,nperm
+ do k=1,ii
+ tabpermchain(k,j)=iequiv(tabperm(k,j),i)
+ enddo
+ enddo
+ else
+ npermc=npermchain
+ npermchain=npermchain*nperm
+ ind=0
+ do k=1,nperm
+ do j=1,npermc
+ ind=ind+1
+ do l=1,ii
+ tabpermchain(l,ind)=tabpermchain(l,j)
+ enddo
+ do l=1,nequiv(i)
+ tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
+ enddo
+ enddo
+ enddo
+ ii=ii+nequiv(i)
+ endif
+ enddo
+ do i=1,npermchain
+ do j=1,nchain
+ itemp(mapchain(j))=tabpermchain(j,i)
+ enddo
+ do j=1,nchain
+ tabpermchain(j,i)=itemp(j)
+ enddo
+ enddo
+ write(iout,*) "Number of chain permutations",npermchain
+ write(iout,*) "Permutations"
+ do i=1,npermchain
+ write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
+ enddo
+ return
+ end
+!c---------------------------------------------------------------------
+ integer function tperm(i,iperm,tabpermchain)
+! implicit none
+! include 'DIMENSIONS'
+ integer i,iperm
+ integer tabpermchain(50,5040)
+ if (i.eq.0) then
+ tperm=0
+ else
+ tperm=tabpermchain(i,iperm)
+ endif
+ return
+ end function
+