From 57fd9101200e85f97d56c97be8c418396b08e05a Mon Sep 17 00:00:00 2001 From: Cezary Czaplewski Date: Sat, 2 May 2020 20:53:04 +0200 Subject: [PATCH] xdrf2pdb Adam's PBC correction --- source/xdrfpdb/src-M/CMakeLists.txt | 6 +- source/xdrfpdb/src-M/COMMON.CHAIN | 12 ++- source/xdrfpdb/src-M/DIMENSIONS | 7 +- source/xdrfpdb/src-M/Makefile | 50 +--------- source/xdrfpdb/src-M/Makefile_gfortran | 54 +++++++++++ source/xdrfpdb/src-M/Makefile_ifort | 54 +++++++++++ source/xdrfpdb/src-M/geomout.F | 3 +- source/xdrfpdb/src-M/seq2chains.f | 56 ++++++++++++ source/xdrfpdb/src-M/xdrf2pdb-m.F | 136 +++++++++++++++++++++++++++- source/xdrfpdb/src-M/xdrf2pdb.F | 156 +++++++++++++++++++++++++++++--- source/xdrfpdb/src-M/xdrf2pdb1.f | 140 ---------------------------- source/xdrfpdb/src-M/xdrf2x1.f | 67 -------------- 12 files changed, 462 insertions(+), 279 deletions(-) mode change 100644 => 120000 source/xdrfpdb/src-M/Makefile create mode 100644 source/xdrfpdb/src-M/Makefile_gfortran create mode 100644 source/xdrfpdb/src-M/Makefile_ifort create mode 100644 source/xdrfpdb/src-M/seq2chains.f delete mode 100644 source/xdrfpdb/src-M/xdrf2pdb1.f delete mode 100644 source/xdrfpdb/src-M/xdrf2x1.f diff --git a/source/xdrfpdb/src-M/CMakeLists.txt b/source/xdrfpdb/src-M/CMakeLists.txt index 2a4ce69..bd9ad24 100644 --- a/source/xdrfpdb/src-M/CMakeLists.txt +++ b/source/xdrfpdb/src-M/CMakeLists.txt @@ -11,8 +11,8 @@ # xdrf2ang : extracts backbone angles from a cx trajectory file # xdrf2pdb-m : converts a selected trajectory of a MREMD run dumpend into a cx file to PDB format # -# xdrf2pdb1 : converts conformation(s) selected from a wham post-processing run into PDB format -# xdrf2x1 : converts conformation(s) selected from a wham post-processing run into raw (x) format. +# removed xdrf2pdb1 : converts conformation(s) selected from a wham post-processing run into PDB format +# removed xdrf2x1 : converts conformation(s) selected from a wham post-processing run into raw (x) format. # # 9/23/2010 A. Liwo # @@ -25,6 +25,7 @@ set(UNRES_XDRF_XDRF2PDB_SRC-M misc.f rescode.f nazwy.f + seq2chains.f ) set(UNRES_XDRF_XDRF2PDB-M_SRC-M @@ -33,6 +34,7 @@ set(UNRES_XDRF_XDRF2PDB-M_SRC-M misc.f rescode.f nazwy.f + seq2chains.f ) diff --git a/source/xdrfpdb/src-M/COMMON.CHAIN b/source/xdrfpdb/src-M/COMMON.CHAIN index 75e8de5..0f42aba 100644 --- a/source/xdrfpdb/src-M/COMMON.CHAIN +++ b/source/xdrfpdb/src-M/COMMON.CHAIN @@ -1,11 +1,19 @@ - integer nres,nsup,nstart_sup,nz_start,nz_end,nres0,nstart_seq + integer nres,nsup,nstart_sup,nz_start,nz_end,nres0,nstart_seq, + & nchain,chain_length,chain_border,iprzes, + & chain_border1,ireschain,tabpermchain,npermchain double precision c,dc,dc_old,d_c_work,xloc,xrot,dc_norm,t,r, & prod,rt,dc_work,cref,crefjlee common /chain/ c(3,maxres2+2),dc(3,0:maxres2),dc_old(3,0:maxres2), & xloc(3,maxres),xrot(3,maxres),dc_norm(3,0:maxres2), - & dc_work(MAXRES6),nres,nres0 + & dc_work(MAXRES6),nres,nres0, + & chain_length(maxchain),npermchain,ireschain(maxres), + & tabpermchain(maxchain,maxperm), + & chain_border(2,maxchain),chain_border1(2,maxchain),nchain common /rotmat/ t(3,3,maxres),r(3,3,maxres),prod(3,3,maxres), & rt(3,3,maxres) common /refstruct/ cref(3,maxres2+2),crefjlee(3,maxres2+2), & nsup,nstart_sup,nstart_seq common /from_zscore/ nz_start,nz_end + double precision boxsize(3) + common /box/ boxsize + diff --git a/source/xdrfpdb/src-M/DIMENSIONS b/source/xdrfpdb/src-M/DIMENSIONS index 43e77a2..3bfc799 100644 --- a/source/xdrfpdb/src-M/DIMENSIONS +++ b/source/xdrfpdb/src-M/DIMENSIONS @@ -15,7 +15,7 @@ C Max. number of coarse-grain processors parameter (max_cg_procs=maxprocs) C Max. number of AA residues integer maxres - parameter (maxres=1500) + parameter (maxres=5000) C Appr. max. number of interaction sites integer maxres2,maxres6,mmaxres6 parameter (maxres2=2*maxres,maxres6=6*maxres) @@ -132,3 +132,8 @@ C Maximum number of conformation stored in cache on each CPU before sending C to master; depends on nstex / ntwx ratio integer max_cache_traj parameter (max_cache_traj=10) +C Max number of symetric chains + integer maxchain + parameter (maxchain=50) + integer maxperm + parameter (maxperm=120) diff --git a/source/xdrfpdb/src-M/Makefile b/source/xdrfpdb/src-M/Makefile deleted file mode 100644 index 72f9dbe..0000000 --- a/source/xdrfpdb/src-M/Makefile +++ /dev/null @@ -1,49 +0,0 @@ -# Set of programs to convert UNRES xdrf format (compressed Cartesian coordinates) to PDF -# or raw-Cartesian format (*.x) or to extract backbone angular coordinates (*.ang) -# The pdb files can be constructed from canonical or MREMD trajectories. -# -# The xdrf library is required -# -# Programs -# -# xdrf2pdb : converts a single cx trajectory file to PDB format -# xdrf2x : converts a single cx trajectory file to raw-coordinate (x) format -# xdrf2ang : extracts backbone angles from a cx trajectory file -# xdrf2pdb-m : converts a selected trajectory of a MREMD run dumpend into a cx file to PDB format -# -# xdrf2pdb1 : converts conformation(s) selected from a wham post-processing run into PDB format -# xdrf2x1 : converts conformation(s) selected from a wham post-processing run into raw (x) format. -# -# 9/23/2010 A. Liwo - -FC=gfortran - -BINDIR = ../../../bin - -#OPT = -fast -OPT = - -FFLAGS = -c ${OPT} - -CPPFLAGS = -DLINUX -DUNRES -DMP -DMPI -DSPLITELE -DPROCOR -DNEWUNRES - -M4 = m4 -M4FILE = underscore.m4 - -.SUFFIXES: .F -.F.o: - ${FC} ${FFLAGS} ${CPPFLAGS} $*.F -.c.o: - ${CC} -c ${CPPFLAGS} $*.c - -xdrf2pdb: xdrf2pdb.o geomout.o misc.o rescode.o nazwy.o xdrf/libxdrf.a - ${FC} -Bstatic -o ${BINDIR}/xdrf2pdb-mult xdrf2pdb.o geomout.o rescode.o misc.o nazwy.o xdrf/libxdrf.a - -xdrf2pdb-m: xdrf2pdb-m.o geomout.o misc.o rescode.o nazwy.o xdrf/libxdrf.a - ${FC} -Bstatic -o ${BINDIR}/xdrf2pdb-m-mult xdrf2pdb-m.o geomout.o rescode.o misc.o nazwy.o xdrf/libxdrf.a - -xdrf2pdb1: xdrf2pdb1.o geomout.o misc.o rescode.o nazwy.o xdrf/libxdrf.a - ${FC} -Bstatic -o ${BINDIR}/xdrf2pdb1-mult xdrf2pdb1.o geomout.o rescode.o misc.o nazwy.o xdrf/libxdrf.a - -clean: - rm -f *.o diff --git a/source/xdrfpdb/src-M/Makefile b/source/xdrfpdb/src-M/Makefile new file mode 120000 index 0000000..696c70e --- /dev/null +++ b/source/xdrfpdb/src-M/Makefile @@ -0,0 +1 @@ +Makefile_ifort \ No newline at end of file diff --git a/source/xdrfpdb/src-M/Makefile_gfortran b/source/xdrfpdb/src-M/Makefile_gfortran new file mode 100644 index 0000000..fe64117 --- /dev/null +++ b/source/xdrfpdb/src-M/Makefile_gfortran @@ -0,0 +1,54 @@ +# Set of programs to convert UNRES xdrf format (compressed Cartesian coordinates) to PDF +# or raw-Cartesian format (*.x) or to extract backbone angular coordinates (*.ang) +# The pdb files can be constructed from canonical or MREMD trajectories. +# +# The xdrf library is required +# +# Programs +# +# xdrf2pdb : converts a single cx trajectory file to PDB format +# xdrf2x : converts a single cx trajectory file to raw-coordinate (x) format +# xdrf2ang : extracts backbone angles from a cx trajectory file +# xdrf2pdb-m : converts a selected trajectory of a MREMD run dumpend into a cx file to PDB format +# +# xdrf2pdb1 : converts conformation(s) selected from a wham post-processing run into PDB format +# xdrf2x1 : converts conformation(s) selected from a wham post-processing run into raw (x) format. +# +# 9/23/2010 A. Liwo + +FC=gfortran + +BINDIR = ../../../bin + +#OPT = -fast +OPT = + +FFLAGS = -c ${OPT} + +CPPFLAGS = -DLINUX -DUNRES -DMP -DMPI -DSPLITELE -DPROCOR -DNEWUNRES + +M4 = m4 +M4FILE = underscore.m4 + +.SUFFIXES: .F +.F.o: + ${FC} ${FFLAGS} ${CPPFLAGS} $*.F +.f.o: + ${FC} ${FFLAGS} ${CPPFLAGS} $*.f +.c.o: + ${CC} -c ${CPPFLAGS} $*.c + +xdrf2pdb: xdrf2pdb.o geomout.o misc.o rescode.o nazwy.o seq2chains.o xdrf/libxdrf.a + ${FC} -Bstatic -o ${BINDIR}/xdrf2pdb-mult xdrf2pdb.o geomout.o rescode.o misc.o nazwy.o seq2chains.o xdrf/libxdrf.a + +xdrf2pdb-m: xdrf2pdb-m.o geomout.o misc.o rescode.o nazwy.o seq2chains.o xdrf/libxdrf.a + ${FC} -Bstatic -o ${BINDIR}/xdrf2pdb-m-mult xdrf2pdb-m.o geomout.o rescode.o misc.o nazwy.o seq2chains.o xdrf/libxdrf.a + +xdrf2x: xdrf2x.o xdrf/libxdrf.a + ${FC} -o xdrf2x xdrf2x-mult.o xdrf/libxdrf.a + +xdrf/libxdrf.a: + cd xdrf && make + +clean: + rm -f *.o diff --git a/source/xdrfpdb/src-M/Makefile_ifort b/source/xdrfpdb/src-M/Makefile_ifort new file mode 100644 index 0000000..cb8d6b1 --- /dev/null +++ b/source/xdrfpdb/src-M/Makefile_ifort @@ -0,0 +1,54 @@ +# Set of programs to convert UNRES xdrf format (compressed Cartesian coordinates) to PDF +# or raw-Cartesian format (*.x) or to extract backbone angular coordinates (*.ang) +# The pdb files can be constructed from canonical or MREMD trajectories. +# +# The xdrf library is required +# +# Programs +# +# xdrf2pdb : converts a single cx trajectory file to PDB format +# xdrf2x : converts a single cx trajectory file to raw-coordinate (x) format +# xdrf2ang : extracts backbone angles from a cx trajectory file +# xdrf2pdb-m : converts a selected trajectory of a MREMD run dumpend into a cx file to PDB format +# +# xdrf2pdb1 : converts conformation(s) selected from a wham post-processing run into PDB format +# xdrf2x1 : converts conformation(s) selected from a wham post-processing run into raw (x) format. +# +# 9/23/2010 A. Liwo + +FC=ifort + +BINDIR = /users2/adam/bin + +#OPT = -fast +OPT = + +FFLAGS = -c ${OPT} + +CPPFLAGS = -DLINUX -DUNRES -DMP -DMPI -DSPLITELE -DPROCOR -DNEWUNRES + +M4 = m4 +M4FILE = underscore.m4 + +.SUFFIXES: .F +.F.o: + ${FC} ${FFLAGS} ${CPPFLAGS} $*.F +.f.o: + ${FC} ${FFLAGS} ${CPPFLAGS} $*.f +.c.o: + ${CC} -c ${CPPFLAGS} $*.c + +xdrf2pdb: xdrf2pdb.o geomout.o misc.o rescode.o nazwy.o seq2chains.o xdrf/libxdrf.a + ${FC} -Bstatic -o ${BINDIR}/xdrf2pdb-mult xdrf2pdb.o geomout.o rescode.o misc.o nazwy.o seq2chains.o xdrf/libxdrf.a + +xdrf2pdb-m: xdrf2pdb-m.o geomout.o misc.o rescode.o nazwy.o seq2chains.o xdrf/libxdrf.a + ${FC} -Bstatic -o ${BINDIR}/xdrf2pdb-m-mult xdrf2pdb-m.o geomout.o rescode.o misc.o nazwy.o seq2chains.o xdrf/libxdrf.a + +xdrf2x: xdrf2x.o xdrf/libxdrf.a + ${FC} -o xdrf2x xdrf2x-mult.o xdrf/libxdrf.a + +xdrf/libxdrf.a: + cd xdrf && make + +clean: + rm -f *.o diff --git a/source/xdrfpdb/src-M/geomout.F b/source/xdrfpdb/src-M/geomout.F index 4b871b8..256bc7a 100644 --- a/source/xdrfpdb/src-M/geomout.F +++ b/source/xdrfpdb/src-M/geomout.F @@ -98,7 +98,7 @@ c endif if (iti.eq.ntyp1) then ichain=ichain+1 ires=0 - write (iunit,'(a)') 'TER' + if (iti_prev.ne.ntyp1) write (iunit,'(a)') 'TER' else ires=ires+1 iatom=iatom+1 @@ -113,6 +113,7 @@ c endif ! & vtot(i+nres) endif endif + iti_prev=iti enddo write (iunit,'(a)') 'TER' do i=nnt,nct-1 diff --git a/source/xdrfpdb/src-M/seq2chains.f b/source/xdrfpdb/src-M/seq2chains.f new file mode 100644 index 0000000..cf38c87 --- /dev/null +++ b/source/xdrfpdb/src-M/seq2chains.f @@ -0,0 +1,56 @@ + subroutine seq2chains(nres,itype,nchain,chain_length,chain_border, + & ireschain) +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' + integer nres,itype(nres),nchain,chain_length(nres), + & chain_border(2,nres),ireschain(nres) + integer ii,ichain,i,j + logical new_chain + ichain=1 + new_chain=.true. + chain_length(ichain)=0 + ii=1 + do while (ii.lt.nres) + if (itype(ii).eq.ntyp1) 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).eq.ntyp1) 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 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 diff --git a/source/xdrfpdb/src-M/xdrf2pdb.F b/source/xdrfpdb/src-M/xdrf2pdb.F index 503d7d7..1286f99 100644 --- a/source/xdrfpdb/src-M/xdrf2pdb.F +++ b/source/xdrfpdb/src-M/xdrf2pdb.F @@ -6,6 +6,7 @@ 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 @@ -17,14 +18,17 @@ 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) @@ -63,26 +67,42 @@ 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 @@ -127,8 +147,6 @@ 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) @@ -149,6 +167,116 @@ c write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize) 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) diff --git a/source/xdrfpdb/src-M/xdrf2pdb1.f b/source/xdrfpdb/src-M/xdrf2pdb1.f deleted file mode 100644 index 7409d7d..0000000 --- a/source/xdrfpdb/src-M/xdrf2pdb1.f +++ /dev/null @@ -1,140 +0,0 @@ - implicit real*8 (a-h,o-z) - include 'DIMENSIONS' - include 'COMMON.CHAIN' - include 'COMMON.INTERACT' - include 'COMMON.SBRIDGE' - real*4 coord(3,1000) - real*4 prec,potE,efree,rmsdev - real*8 etot - character*80 arg,seqfile,pdbfile - character*3 sequenc(maxres) - character*50 tytul - character*8 onethree,cfreq - character*8 ucase - external ucase - logical oneletter - integer rescode - external rescode - logical ilen,iblnk - - ifreq=1 - if (iargc().lt.4) then - print '(a)', - & "Usage: xdrf2pdb one/three seqfile cxfile conf [pdbfile]" - 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 - print *,"nres",nres - do i=1,nres - print *,i," ",sequenc(i) - enddo - print * - i=0 - do while (.not.iblnk(sequenc(i+1)(1:1))) - print *,i+1," ",sequenc(i+1)," ",sequenc(i+1)(1:1) - i=i+1 - enddo - nres=i - print *,"nres",nres - do i=1,nres - itype(i)=rescode(i,sequenc(i),0) - enddo - print *,(itype(i),i=1,nres) - endif - call getarg(3,arg) - iext = index(arg,'.cx') - 1 - if (iext.lt.0) then - print *,"Error - not a cx file" - stop - endif - print *,"arg ",arg - call getarg(4,cfreq) - read (cfreq,*) iconf - print *,"iconf",iconf - if (iargc().gt.4) then - call getarg(5,pdbfile) - else - pdbfile=arg(:iext)//'.pdb' - endif - open(9,file=pdbfile) - nnt = 1 - if (itype(1).eq.21) nnt = 2 - nct=nres - if (itype(nres).eq.21) nct = nres-1 -c if (nct.eq.nres-1) nres=nres-1 -c if (nnt.eq.2) nres=nres-1 - - print *,"nres",nres," nnt",nnt," nct",nct - - call xdrfopen(ixdrf,arg, "r", iret) - print *,"iret",iret - kk = 0 - do while(.true.) - prec=10000.0 - isize=0 - call xdrf3dfcoord(ixdrf, coord, isize, prec, iret) - call xdrfint(ixdrf, nss, iret) - do j=1,nss - call xdrfint(ixdrf, ihpb(j), iret) - call xdrfint(ixdrf, jhpb(j), iret) - enddo - call xdrffloat(ixdrf, potE, iret) - if(iret.eq.0) exit - kk = kk + 1 - call xdrffloat(ixdrf, efree, iret) - call xdrffloat(ixdrf, rmsdev, iret) - call xdrfint(ixdrf, iscor, iret) - if (kk.eq.iconf) then - print *,"pote",pote," efree",efree," rmsdev",rmsdev - - print *,"isize",isize - - 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 - etot=potE - write (tytul,'(a,i6)') "Structure",kk - call pdbout(etot,tytul,9) - stop - endif - enddo - - end diff --git a/source/xdrfpdb/src-M/xdrf2x1.f b/source/xdrfpdb/src-M/xdrf2x1.f deleted file mode 100644 index 53f8fb0..0000000 --- a/source/xdrfpdb/src-M/xdrf2x1.f +++ /dev/null @@ -1,67 +0,0 @@ - implicit real*8 (a-h,o-z) - integer ihpb(100),jhpb(100) - real*4 coord(3,1000) - real*4 prec,potE,efree,rmsdev,qfrag(100) - real*8 etot - character*80 arg,xfile - character*8 ucase,cfreq - external ucase - logical oneletter - integer rescode - external rescode - - ifreq=1 - if (iargc().lt.2) then - print '(a)', - & "Usage: xdrf2x1 cxfile conf [pdbfile]" - stop - endif - call getarg(1,arg) - iext = index(arg,'.cx') - 1 - if (iext.lt.0) then - print *,"Error - not a cx file" - stop - endif - print *,"arg ",arg - call getarg(2,cfreq) - read (cfreq,*) iconf - print *,"iconf",iconf - if (iargc().gt.2) then - call getarg(3,pdbfile) - else - xfile=arg(:iext)//'.x' - endif - open(9,file=xfile) - - call xdrfopen(ixdrf,arg, "r", iret) - print *,"iret",iret - kk = 0 - do while(.true.) - prec=10000.0 - isize=0 - call xdrf3dfcoord(ixdrf, coord, isize, prec, iret) - call xdrfint(ixdrf, nss, iret) - do j=1,nss - call xdrfint(ixdrf, ihpb(j), iret) - call xdrfint(ixdrf, jhpb(j), iret) - enddo - call xdrffloat(ixdrf, potE, iret) - if(iret.eq.0) exit - kk = kk + 1 - call xdrffloat(ixdrf, efree, iret) - call xdrffloat(ixdrf, rmsdev, iret) - call xdrfint(ixdrf, iscor, iret) - if (kk.eq.iconf) then - print *,"pote",pote," efree",efree," rmsdev",rmsdev - - print *,"isize",isize - - write (9,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath - write (9,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss) - write (9,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag) - write (9,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize) - stop - endif - enddo - - end -- 1.7.9.5