NEWCORR5D working with 15k, work and iwork in random_vel might need testing
[unres4.git] / source / unres / io.F90
index 77bd83a..4750e57 100644 (file)
       use MPI_data
       use compare, only:seq_comp,contact
       use control
+      implicit none
 !      implicit real*8 (a-h,o-z)
 !      include 'DIMENSIONS'
 #ifdef MPI
 
 !        print *,'Finished reading pdb data'
         if(me.eq.king.or..not.out1file) &
-         write (iout,'(a,i3,a,i3)')'nsup=',nsup,&
+         write (iout,'(a,i7,a,i7)')'nsup=',nsup,&
          ' nstart_sup=',nstart_sup !,"ergwergewrgae"
 !el        if(.not.allocated(itype_pdb)) 
         allocate(itype_pdb(nres))
 !       print *,"tu dochodze"
 ! znamy nres oraz nss można zaalokowac potrzebne tablice
       call alloc_geo_arrays
+      print *,"after  alloc_geo_arrays"
       call alloc_ener_arrays
+      print *,"after  alloc_ener_arrays"
 !--------------------------------
 ! 8/13/98 Set limits to generating the dihedral angles
       do i=1,nres
       endif
 #endif
       nct=nres
+      allocate(ireschain(nres))
+      ireschain=0
+      write(iout,*),"before seq2chains",ireschain
+      call seq2chains
+      write(iout,*) "after seq2chains",nchain
+      allocate ( chain_border1(2,nchain))
+      chain_border1(1,1)=1
+      chain_border1(2,1)=chain_border(2,1)+1
+      do i=2,nchain-1
+        chain_border1(1,i)=chain_border(1,i)-1
+        chain_border1(2,i)=chain_border(2,i)+1
+      enddo
+      if (nchain.gt.1) chain_border1(1,nchain)=chain_border(1,nchain)-1
+      chain_border1(2,nchain)=nres
+      write(iout,*) "nres",nres," nchain",nchain
+      do i=1,nchain
+        write(iout,*)"chain",i,chain_length(i),chain_border(1,i),&
+         chain_border(2,i),chain_border1(1,i),chain_border1(2,i)
+      enddo
+      allocate(tabpermchain(50,5040))
+      call chain_symmetry(npermchain,tabpermchain)
       print *,'NNT=',NNT,' NCT=',NCT
       if (itype(1,molnum(1)).eq.ntyp1_molec(molnum(1))) nnt=2
       if (itype(nres,molnum(nres)).eq.ntyp1_molec(molnum(nres))) nct=nct-1
       if (pdbref) then
         if(me.eq.king.or..not.out1file) &
-         write (iout,'(a,i3)') 'nsup=',nsup
+         write (iout,'(a,i7)') 'nsup=',nsup
         nstart_seq=nnt
         if (nsup.le.(nct-nnt+1)) then
           do i=0,nct-nnt+1-nsup
       return
    10 stop "Error in fragment file"
       end subroutine read_klapaucjusz
+
+!-----------------------------------------------------------
+      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
+
 !-----------------------------------------------------------------------------
       end module io