X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fxdrfpdb%2Fsrc-M%2Fxdrf2pdb-m.F;fp=source%2Fxdrfpdb%2Fsrc-M%2Fxdrf2pdb-m.F;h=e3825b08917048f017885d7ea0cbafd2f33d39a4;hb=0f8124c6b34a8a30b39839f0e4bf8c6c20dcbb4e;hp=d2d7d9c4c5c09bc6a0abc8dcb0413ff4cd382d36;hpb=5184115e0e952c87579f5dde1a6ad30d28eb0ea1;p=unres.git diff --git a/source/xdrfpdb/src-M/xdrf2pdb-m.F b/source/xdrfpdb/src-M/xdrf2pdb-m.F index d2d7d9c..e3825b0 100644 --- a/source/xdrfpdb/src-M/xdrf2pdb-m.F +++ b/source/xdrfpdb/src-M/xdrf2pdb-m.F @@ -16,12 +16,15 @@ logical oneletter,iblnk integer rescode external rescode + integer lshift(maxchain),lmin(maxchain) + double precision shift_coord(3,maxchain),cm(3,maxchain),ccm(3), + & ccm_prev(3) ifreq=1 - if (iargc().lt.5) then + if (iargc().lt.8) then print '(2a)', - & "Usage: xdrf2pdb-m one/three seqfile cxfile ntraj itraj", - & " [pdbfile] [freq]" + & "Usage: xdrf2pdb-m one/three seqfile boxx boxy boxz cxfile ", + & "ntraj itraj [pdbfile] [freq]" stop endif call getarg(1,onethree) @@ -62,6 +65,22 @@ print *,nres print '(20(a3,1x))',(sequenc(i),i=1,nres) 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(3,arg) iext = index(arg,'.cx') - 1 if (iext.lt.0) then @@ -137,6 +156,117 @@ c write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize) c(j,i+nres)=coord(j,ii+nres) enddo enddo +c Bring all chains to same box and to a common origin + 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,'(a1,i6,a8,f6.1,a6,f10.1)') "#",kk," t_bath ", & t_bath, " time ",time