X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;ds=sidebyside;f=source%2Funres%2Fio.F90;h=4750e579ee861416b994ecdb6493941cbf6c9971;hb=91b9f78d94b96277537615722323ebe03cc0a014;hp=92536385731340981cc0ba3073756548d1aad793;hpb=ae2f9918e1e5f88a78cc32b1d3a7ff90dede1207;p=unres4.git diff --git a/source/unres/io.F90 b/source/unres/io.F90 index 9253638..4750e57 100644 --- a/source/unres/io.F90 +++ b/source/unres/io.F90 @@ -712,6 +712,7 @@ use MPI_data use compare, only:seq_comp,contact use control + implicit none ! implicit real*8 (a-h,o-z) ! include 'DIMENSIONS' #ifdef MPI @@ -1014,7 +1015,7 @@ ! 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)) @@ -1154,6 +1155,7 @@ endif enddo do i=2,nres-1 + mnum=molnum(i) if (molnum(i).eq.1) then ! print *, "molnum",molnum(i),itype(i,molnum(i)),nres,i vbld(i+nres)=dsc(iabs(itype(i,molnum(i)))) @@ -1237,7 +1239,9 @@ ! 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 @@ -1368,12 +1372,33 @@ 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 @@ -2629,5 +2654,205 @@ 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