Adam's changes
[unres.git] / source / wham / src-M-SAXS / 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