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