real*4 coord(3,2*maxres)
real*4 prec,time,potE,uconst,t_bath,qfrag(100)
real*8 etot
+ real*8 cm(3,maxchain),shift_coord(3,maxchain),ccm(3),ccm_prev(3)
character*80 arg,seqfile,pdbfile
character*3 sequenc(maxres)
character*50 tytul
external rescode
logical iblnk
external iblnk
+ logical previous,change
+ integer lmin(maxchain),lshift(maxchain)
+ previous=.false.
ifreq=1
is=1
ie=1000000000
- if (iargc().lt.3) then
- print '(2a)',
- & "Usage: xdrf2pdb one/three seqfile cxfile [freq]",
- & " [start] [end] [pdbfile]"
+ if (iargc().lt.6) then
+ print '(2a)',
+ & "Usage: xdrf2pdb one/three seqfile boxx boxy boxz cxfile [freq]",
+ &" [start] [end] [pdbfile]"
stop
endif
call getarg(1,onethree)
itype(i)=rescode(i,sequenc(i),0)
enddo
endif
+ call seq2chains(nres,itype,nchain,chain_length,chain_border,
+ & ireschain)
+ print *,"nchain",nchain
+ print *,"chain_border, chain_length"
+ nrestot=0
+ do i=1,nchain
+ nrestot=nrestot+chain_length(i)
+ print *,chain_border(1,i),chain_border(2,i),chain_length(i)
+ enddo
+ print *,"nrestot",nrestot
call getarg(3,arg)
+ read(arg,*) boxsize(1)
+ call getarg(4,arg)
+ read(arg,*) boxsize(2)
+ call getarg(5,arg)
+ read(arg,*) boxsize(3)
+ call getarg(6,arg)
iext = index(arg,'.cx') - 1
if (iext.lt.0) then
print *,"Error - not a cx file"
stop
endif
- if (iargc().gt.3) then
- call getarg(4,cfreq)
+ if (iargc().gt.6) then
+ call getarg(7,cfreq)
read (cfreq,*) ifreq
endif
- if (iargc().gt.4) then
- call getarg(5,cfreq)
+ if (iargc().gt.7) then
+ call getarg(8,cfreq)
read (cfreq,*) is
endif
- if (iargc().gt.5) then
- call getarg(6,cfreq)
+ if (iargc().gt.8) then
+ call getarg(9,cfreq)
read (cfreq,*) ie
endif
- if (iargc().gt.6) then
- call getarg(7,pdbfile)
+ if (iargc().gt.9) then
+ call getarg(10,pdbfile)
else
pdbfile=arg(:iext)//'.pdb'
endif
isize=0
call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
-
-
c write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
c write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
c write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
c(j,i+nres)=coord(j,ii+nres)
enddo
enddo
+ do ichain=1,nchain
+c print *,"ichain",ichain
+ do j=1,3
+ cm(j,ichain)=0.0d0
+ enddo
+ do i=chain_border(1,ichain),chain_border(2,ichain)
+c print *,i,c(:,i)
+ do j=1,3
+ cm(j,ichain)=cm(j,ichain)+c(j,i)
+ enddo
+ enddo
+ do j=1,3
+ cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
+ enddo
+ enddo
+c write (*,*) "CM"
+c do i=1,nchain
+c write (*,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
+c enddo
+ do j=1,3
+ lmin=0
+ sumodl_min=1.0d10
+ do i=0,3**nchain-1
+ ii=i
+ do ichain=1,nchain
+ lshift(ichain)=mod(ii,3)-1
+ ii=ii/3
+ enddo
+c print '(10i3)',(lshift(ichain),ichain=1,nchain)
+ sumodl=0.0d0
+ do ichain=1,nchain
+ do jchain=1,ichain-1
+ sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
+ & -cm(j,jchain)-lshift(jchain)*boxsize(j))
+ enddo
+ enddo
+c print *,"sumodl",sumodl," sumodl_min",sumodl_min
+ if (sumodl.lt.sumodl_min) then
+ sumodl_min=sumodl
+ lmin=lshift
+ endif
+ enddo
+ do ichain=1,nchain
+ cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
+ shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
+ enddo
+ enddo
+ do ichain=1,nchain
+ do j=1,3
+ ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
+ enddo
+ enddo
+ do j=1,3
+ ccm(j)=ccm(j)/nrestot
+ enddo
+c write (*,'(a,3f10.5)') "ccm ",(ccm(j),j=1,3)
+c write (*,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
+c write (*,'(a,3f10.5)') "diff ",(ccm(j)-ccm_prev(j),j=1,3)
+c write (*,*) "shift_coord"
+c do i=1,nchain
+c write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+c enddo
+#ifdef PREVIOUS
+ if (previous) then
+ do j=1,3
+ change=.true.
+ do while (change)
+ if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
+ do ichain=1,nchain
+ shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
+ enddo
+ ccm(j)=ccm(j)-boxsize(j)
+ else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
+ do ichain=1,nchain
+ shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
+ enddo
+ ccm(j)=ccm(j)+boxsize(j)
+ else
+ change=.false.
+ endif
+ enddo
+ enddo
+ do j=1,3
+ ccm_prev(j)=ccm(j)
+ enddo
+ else
+ previous=.true.
+ do j=1,3
+ ccm_prev(j)=ccm(j)
+ enddo
+ endif
+#else
+ do ichain=1,nchain
+ do j=1,3
+ shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
+ enddo
+ enddo
+#endif
+c write (*,*) "shift_coord"
+c do i=1,nchain
+c write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
+c enddo
+ do ichain=1,nchain
+ do i=chain_border(1,ichain),chain_border(2,ichain)
+ do j=1,3
+ c(j,i)=c(j,i)+shift_coord(j,ichain)
+ c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
+ enddo
+ enddo
+ enddo
etot=potE
write (tytul,'(a,i6)') "Structure",kk
call pdbout(etot,tytul,9)
projects / unres.git / blobdiff