implicit real*8 (a-h,o-z) include 'DIMENSIONS' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.SBRIDGE' real*4 coord(3,5000) real*4 prec,time,potE,uconst,t_bath,qfrag(100) real*8 etot character*80 arg,seqfile,pdbfile character*3 sequenc(maxres) character*50 tytul character*8 onethree,cfreq,cntraj,citraj character*3 licz character*8 ucase external ucase 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.8) then print '(2a)', & "Usage: xdrf2pdb-m one/three seqfile boxx boxy boxz cxfile ", & "ntraj itraj [pdbfile] [freq]" stop endif call getarg(1,onethree) onethree = ucase(onethree) if (onethree.eq.'ONE') then oneletter = .true. else if (onethree.eq.'THREE') then oneletter = .false. else print *,"ONE or THREE must be specified" endif call getarg(2,seqfile) open (1,file=seqfile,status='old') if (oneletter) then read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres) 10 continue nres=i i=0 do while (.not.iblnk(sequenc(i+1)(1:1))) i=i+1 enddo nres=i do i=1,nres itype(i)=rescode(i,sequenc(i),1) enddo else read(1,'(20(a3,1x))',end=11,err=11) (sequenc(i),i=1,maxres) 11 continue nres=i i=0 do while (.not.iblnk(sequenc(i+1)(1:1))) i=i+1 enddo nres=i do i=1,nres itype(i)=rescode(i,sequenc(i),0) enddo 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(6,arg) iext = index(arg,'.cx') - 1 if (iext.lt.0) then print *,"Error - not a cx file" print *,arg stop endif call getarg(7,cntraj) read (cntraj,*) ntraj call getarg(8,citraj) read (citraj,*) itraj if (iargc().gt.8) then call getarg(9,pdbfile) else write(licz,'(bz,i3.3)') itraj pdbfile=arg(:iext)//'_'//licz//'.pdb' endif if (iargc().gt.9) then call getarg(10,cfreq) read (cfreq,*) ifreq endif c print *,"ifreq",ifreq," ntraj",ntraj," itraj",itraj open(9,file=pdbfile) nnt = 1 if (itype(1).eq.ntyp1) nnt = 2 nct=nres if (itype(nres).eq.ntyp1) nct = nres-1 print *,"nnt",nnt," nct",nct call xdrfopen(ixdrf,arg, "r", iret) kk = 0 do while(.true.) call xdrffloat(ixdrf, time, iret) if(iret.eq.0) exit kk = kk + 1 call xdrffloat(ixdrf, potE, iret) call xdrffloat(ixdrf, uconst, iret) c#ifdef NEWUNRES c call xdrffloat(ixdrf, uconst_back, iret) c#endif call xdrffloat(ixdrf, t_bath, iret) call xdrfint(ixdrf, nss, iret) do j=1,nss call xdrfint(ixdrf, ihpb(j), iret) call xdrfint(ixdrf, jhpb(j), iret) enddo call xdrfint(ixdrf, nfrag, iret) call xdrfint(ixdrf, iset, iret) do i=1,nfrag call xdrffloat(ixdrf, qfrag(i), iret) enddo prec=10000.0 c print *,'nss=',nss,'nfrag=',nfrag,'iset=',iset 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 write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize) if (mod(kk/ntraj,ifreq).eq.0 .and. mod(kk,ntraj).eq.itraj) then if (isize .ne. nres+nct-nnt+1) then print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1 endif do i=1,nres do j=1,3 c(j,i)=coord(j,i) enddo enddo ii = 0 do i=nnt,nct ii = ii + 1 do j=1,3 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 call pdbout(etot,tytul,9) endif enddo end