source/xdrfpdb/src-M/xdrf2pdb.F

   1       implicit real*8 (a-h,o-z)
   2       include 'DIMENSIONS'
   3       include 'COMMON.CHAIN'
   4       include 'COMMON.INTERACT'
   5       include 'COMMON.SBRIDGE'
   6       real*4 coord(3,2*maxres)
   7       real*4 prec,time,potE,uconst,t_bath,qfrag(100)
   8       real*8 etot
   9       real*8 cm(3,maxchain),shift_coord(3,maxchain),ccm(3),ccm_prev(3)
  10       character*80 arg,seqfile,pdbfile
  11       character*3 sequenc(maxres)
  12       character*50 tytul
  13       character*8 onethree,cfreq
  14       character*8 ucase
  15       external ucase
  16       logical oneletter
  17       integer rescode
  18       external rescode
  19       logical iblnk
  20       external iblnk
  21       logical previous,change
  22       integer lmin(maxchain),lshift(maxchain)
  23
  24       previous=.false.
  25       ifreq=1
  26       is=1
  27       ie=1000000000
  28       if (iargc().lt.6) then
  29        print '(2a)',
  30      & "Usage: xdrf2pdb one/three seqfile boxx boxy boxz cxfile [freq]",
  31      &" [start] [end] [pdbfile]"
  32         stop
  33       endif
  34       call getarg(1,onethree)
  35       onethree = ucase(onethree)
  36       if (onethree.eq.'ONE') then
  37         oneletter = .true.
  38       else if (onethree.eq.'THREE') then
  39         oneletter = .false.
  40       else
  41         print *,"ONE or THREE must be specified"
  42       endif
  43       call getarg(2,seqfile)
  44       open (1,file=seqfile,status='old')
  45       if (oneletter) then
  46         read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
  47    10   continue
  48         nres=i
  49         i=0
  50         do while (.not.iblnk(sequenc(i+1)(1:1)))
  51           i=i+1
  52         enddo
  53         nres=i
  54         do i=1,nres
  55           itype(i)=rescode(i,sequenc(i),1)
  56         enddo
  57       else
  58         read(1,'(20(a3,1x))',end=11,err=11) (sequenc(i),i=1,maxres)
  59    11   continue
  60         nres=i
  61         i=0
  62         do while (.not.iblnk(sequenc(i+1)(1:1)))
  63           i=i+1
  64         enddo
  65         nres=i
  66         do i=1,nres
  67           itype(i)=rescode(i,sequenc(i),0)
  68         enddo
  69       endif
  70       call seq2chains(nres,itype,nchain,chain_length,chain_border,
  71      &  ireschain)
  72       print *,"nchain",nchain
  73       print *,"chain_border, chain_length"
  74       nrestot=0
  75       do i=1,nchain
  76         nrestot=nrestot+chain_length(i)
  77         print *,chain_border(1,i),chain_border(2,i),chain_length(i)
  78       enddo
  79       print *,"nrestot",nrestot
  80       call getarg(3,arg)
  81       read(arg,*) boxsize(1)
  82       call getarg(4,arg)
  83       read(arg,*) boxsize(2)
  84       call getarg(5,arg)
  85       read(arg,*) boxsize(3)
  86       call getarg(6,arg)
  87       iext = index(arg,'.cx') - 1
  88       if (iext.lt.0) then
  89         print *,"Error - not a cx file"
  90         stop
  91       endif
  92       if (iargc().gt.6) then
  93         call getarg(7,cfreq)
  94         read (cfreq,*) ifreq
  95       endif
  96       if (iargc().gt.7) then
  97         call getarg(8,cfreq)
  98         read (cfreq,*) is
  99       endif
 100       if (iargc().gt.8) then
 101         call getarg(9,cfreq)
 102         read (cfreq,*) ie
 103       endif
 104       if (iargc().gt.9) then
 105         call getarg(10,pdbfile)
 106       else
 107         pdbfile=arg(:iext)//'.pdb'
 108       endif
 109       open(9,file=pdbfile)
 110       nnt = 1
 111       if (itype(1).eq.ntyp1) nnt = 2
 112       nct=nres
 113       if (itype(nres).eq.ntyp1) nct = nres-1
 114       print *,"nnt",nnt," nct",nct
 115       print *,"file",arg
 116       call xdrfopen(ixdrf,arg, "r", iret)
 117       if (iret.eq.0) stop
 118       print *,"iret",iret
 119       print *,"is",is," ie",ie
 120       kk = 0
 121       do while(is.eq.0 .or. kk.lt.ie)
 122        call xdrffloat(ixdrf, time, iret)
 123        print *,"time",time," iret",iret
 124        if(iret.eq.0) exit
 125        kk = kk + 1
 126        call xdrffloat(ixdrf, potE, iret)
 127        call xdrffloat(ixdrf, uconst, iret)
 128        print *,"potE",potE," uconst",uconst
 129 #ifdef NEWUNRES
 130        call xdrffloat(ixdrf, uconst_back, iret)
 131 #endif
 132        print *,"uconst_back",uconst_back
 133        call xdrffloat(ixdrf, t_bath, iret)
 134        print *,"t_bath",t_bath
 135        call xdrfint(ixdrf, nss, iret)
 136        do j=1,nss
 137         call xdrfint(ixdrf, ihpb(j), iret)
 138         call xdrfint(ixdrf, jhpb(j), iret)
 139        enddo
 140        print *,"nss",nss
 141        call xdrfint(ixdrf, nfrag, iret)
 142        do i=1,nfrag
 143         call xdrffloat(ixdrf, qfrag(i), iret)
 144        enddo
 145        print *,"nfrag",nfrag
 146        prec=10000.0
 147
 148        isize=0
 149        call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
 150 c       write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 151 c       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
 152 c       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
 153 c       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
 154        if (kk.ge.is .and. mod(kk,ifreq).eq.0) then
 155          if (isize .ne. nres+nct-nnt+1) then
 156            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
 157          endif
 158          do i=1,nres
 159            do j=1,3
 160              c(j,i)=coord(j,i)
 161            enddo
 162          enddo
 163          ii = 0
 164          do i=nnt,nct
 165            ii = ii + 1
 166            do j=1,3
 167              c(j,i+nres)=coord(j,ii+nres)
 168            enddo
 169          enddo
 170          do ichain=1,nchain
 171 c           print *,"ichain",ichain
 172            do j=1,3
 173              cm(j,ichain)=0.0d0
 174            enddo
 175            do i=chain_border(1,ichain),chain_border(2,ichain)
 176 c             print *,i,c(:,i)
 177              do j=1,3
 178                cm(j,ichain)=cm(j,ichain)+c(j,i)
 179              enddo
 180            enddo
 181            do j=1,3
 182              cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
 183            enddo
 184          enddo
 185 c         write (*,*) "CM"
 186 c         do i=1,nchain
 187 c           write (*,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
 188 c         enddo
 189          do j=1,3
 190            lmin=0
 191            sumodl_min=1.0d10
 192            do i=0,3**nchain-1
 193              ii=i
 194              do ichain=1,nchain
 195                lshift(ichain)=mod(ii,3)-1
 196                ii=ii/3
 197              enddo
 198 c             print '(10i3)',(lshift(ichain),ichain=1,nchain)
 199              sumodl=0.0d0
 200              do ichain=1,nchain
 201                do jchain=1,ichain-1
 202                sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
 203      &             -cm(j,jchain)-lshift(jchain)*boxsize(j))
 204                enddo
 205              enddo
 206 c             print *,"sumodl",sumodl," sumodl_min",sumodl_min
 207              if (sumodl.lt.sumodl_min) then
 208                sumodl_min=sumodl
 209                lmin=lshift
 210              endif
 211            enddo
 212            do ichain=1,nchain
 213              cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
 214              shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
 215            enddo
 216          enddo
 217          do ichain=1,nchain
 218            do j=1,3
 219              ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
 220            enddo
 221          enddo
 222          do j=1,3
 223            ccm(j)=ccm(j)/nrestot
 224          enddo
 225 c         write (*,'(a,3f10.5)') "ccm     ",(ccm(j),j=1,3)
 226 c         write (*,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
 227 c         write (*,'(a,3f10.5)') "diff    ",(ccm(j)-ccm_prev(j),j=1,3)
 228 c         write (*,*) "shift_coord"
 229 c         do i=1,nchain
 230 c           write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
 231 c         enddo
 232 #ifdef PREVIOUS
 233          if (previous) then
 234            do j=1,3
 235              change=.true.
 236              do while (change)
 237              if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
 238                do ichain=1,nchain
 239                  shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
 240                enddo
 241                ccm(j)=ccm(j)-boxsize(j)
 242              else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
 243                do ichain=1,nchain
 244                  shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
 245                enddo
 246                ccm(j)=ccm(j)+boxsize(j)
 247              else
 248                change=.false.
 249              endif
 250            enddo
 251            enddo
 252            do j=1,3
 253              ccm_prev(j)=ccm(j)
 254            enddo
 255          else
 256            previous=.true.
 257            do j=1,3
 258              ccm_prev(j)=ccm(j)
 259            enddo
 260          endif
 261 #else
 262          do ichain=1,nchain
 263            do j=1,3
 264              shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
 265            enddo
 266          enddo
 267 #endif
 268 c         write (*,*) "shift_coord"
 269 c         do i=1,nchain
 270 c           write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
 271 c         enddo
 272          do ichain=1,nchain
 273            do i=chain_border(1,ichain),chain_border(2,ichain)
 274              do j=1,3
 275                c(j,i)=c(j,i)+shift_coord(j,ichain)
 276                c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
 277              enddo
 278            enddo
 279          enddo
 280          etot=potE
 281          write (tytul,'(a,i6)') "Structure",kk
 282          call pdbout(etot,tytul,9)
 283        endif
 284       enddo
 285
 286       end