source/wham/src-HCD-5D/chain_symmetry.F

   1       subroutine chain_symmetry(nchain,nres,itype,chain_border,
   2      &    chain_length,npermchain,tabpermchain)
   3 c
   4 c Determine chain symmetry. nperm is the number of permutations and
   5 c tabperchain contains the allowed permutations of the chains.
   6 c
   7       implicit none
   8       include "DIMENSIONS"
   9       include "COMMON.IOUNITS"
  10       integer nchain,nres,itype(nres),chain_border(2,maxchain),
  11      &  chain_length(nchain),itemp(maxchain),
  12      &  npermchain,tabpermchain(maxchain,maxperm),
  13      &  tabperm(maxchain,maxperm),mapchain(maxchain),
  14      &  iequiv(maxchain,maxres),iflag(maxres)
  15       integer i,j,k,l,ii,nchain_group,nequiv(maxchain),iieq,
  16      &  nperm,npermc,ind
  17       if (nchain.eq.1) then
  18         npermchain=1
  19         tabpermchain(1,1)=1
  20 c        print*,"npermchain",npermchain," tabpermchain",tabpermchain(1,1)
  21         return
  22       endif
  23 c
  24 c Look for equivalent chains
  25 c
  26 #ifdef DEBUG
  27       write(iout,*) "nchain",nchain
  28       do i=1,nchain
  29         write(iout,*) "chain",i," from",chain_border(1,i),
  30      &      " to",chain_border(2,i)
  31         write(iout,*)
  32      &   "sequence ",(itype(j),j=chain_border(1,i),chain_border(2,i))
  33       enddo
  34 #endif
  35       do i=1,nchain
  36         iflag(i)=0
  37       enddo
  38       nchain_group=0
  39       do i=1,nchain
  40         if (iflag(i).gt.0) cycle
  41         iflag(i)=1
  42         nchain_group=nchain_group+1
  43         iieq=1
  44         iequiv(iieq,nchain_group)=i
  45         do j=i+1,nchain
  46           if (iflag(j).gt.0.or.chain_length(i).ne.chain_length(j)) cycle
  47 c          k=0
  48 c          do while(k.lt.chain_length(i) .and.
  49 c     &     itype(chain_border(1,i)+k).eq.itype(chain_border(1,j)+k))
  50           do k=0,chain_length(i)-1
  51 c            k=k+1
  52             if (itype(chain_border(1,i)+k).ne.
  53      &          itype(chain_border(1,j)+k)) exit
  54           enddo
  55           if (k.lt.chain_length(i)) cycle
  56           iflag(j)=1
  57           iieq=iieq+1
  58           iequiv(iieq,nchain_group)=j
  59         enddo
  60         nequiv(nchain_group)=iieq
  61       enddo
  62       write(iout,*) "Number of equivalent chain groups:",nchain_group
  63       write(iout,*) "Equivalent chain groups"
  64       do i=1,nchain_group
  65         write(iout,*) "group",i," #members",nequiv(i)," chains",
  66      &      (iequiv(j,i),j=1,nequiv(i))
  67       enddo
  68       ind=0
  69       do i=1,nchain_group
  70         do j=1,nequiv(i)
  71           ind=ind+1
  72           mapchain(ind)=iequiv(j,i)
  73         enddo
  74       enddo
  75       write (iout,*) "mapchain"
  76       do i=1,nchain
  77         write (iout,*) i,mapchain(i)
  78       enddo
  79       ii=0
  80       do i=1,nchain_group
  81         call permut(nequiv(i),nperm,tabperm)
  82         if (ii.eq.0) then
  83           ii=nequiv(i)
  84           npermchain=nperm
  85           do j=1,nperm
  86             do k=1,ii
  87               tabpermchain(k,j)=iequiv(tabperm(k,j),i)
  88             enddo
  89           enddo
  90         else
  91           npermc=npermchain
  92           npermchain=npermchain*nperm
  93           ind=0
  94           do k=1,nperm
  95             do j=1,npermc
  96               ind=ind+1
  97               do l=1,ii
  98                 tabpermchain(l,ind)=tabpermchain(l,j)
  99               enddo
 100               do l=1,nequiv(i)
 101                 tabpermchain(ii+l,ind)=iequiv(tabperm(l,k),i)
 102               enddo
 103             enddo
 104           enddo
 105           ii=ii+nequiv(i)
 106         endif
 107       enddo
 108       do i=1,npermchain
 109         do j=1,nchain
 110           itemp(mapchain(j))=tabpermchain(j,i)
 111         enddo
 112         do j=1,nchain
 113           tabpermchain(j,i)=itemp(j)
 114         enddo
 115       enddo
 116       write(iout,*) "Number of chain permutations",npermchain
 117       write(iout,*) "Permutations"
 118       do i=1,npermchain
 119         write(iout,'(20i4)') (tabpermchain(j,i),j=1,nchain)
 120       enddo
 121       return
 122       end
 123 c---------------------------------------------------------------------
 124       integer function tperm(i,iperm,tabpermchain)
 125       implicit none
 126       include 'DIMENSIONS'
 127       integer i,iperm
 128       integer tabpermchain(maxchain,maxperm)
 129       if (i.eq.0) then
 130         tperm=0
 131       else
 132         tperm=tabpermchain(i,iperm)
 133       endif
 134       return
 135       end