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(6,arg)
  85       iext = index(arg,'.cx') - 1
  86       if (iext.lt.0) then
  87         print *,"Error - not a cx file"
  88         print *,arg
  89         stop
  90       endif
  91       call getarg(7,cntraj)
  92       read (cntraj,*) ntraj
  93       call getarg(8,citraj)
  94       read (citraj,*) itraj
  95       if (iargc().gt.8) then
  96         call getarg(9,pdbfile)
  97       else
  98         write(licz,'(bz,i3.3)') itraj
  99         pdbfile=arg(:iext)//'_'//licz//'.pdb'
 100       endif
 101       if (iargc().gt.9) then
 102         call getarg(10,cfreq)
 103         read (cfreq,*) ifreq
 104       endif
 105 c      print *,"ifreq",ifreq," ntraj",ntraj," itraj",itraj
 106       open(9,file=pdbfile)
 107       nnt = 1
 108       if (itype(1).eq.ntyp1) nnt = 2
 109       nct=nres
 110       if (itype(nres).eq.ntyp1) nct = nres-1
 111       print *,"nnt",nnt," nct",nct
 112       call xdrfopen(ixdrf,arg, "r", iret)
 113       kk = 0
 114       do while(.true.)
 115        call xdrffloat(ixdrf, time, iret)
 116        if(iret.eq.0) exit
 117        kk = kk + 1
 118        call xdrffloat(ixdrf, potE, iret)
 119        call xdrffloat(ixdrf, uconst, iret)
 120 c#ifdef NEWUNRES
 121 c       call xdrffloat(ixdrf, uconst_back, iret)
 122 c#endif
 123        call xdrffloat(ixdrf, t_bath, iret)
 124        call xdrfint(ixdrf, nss, iret)
 125        do j=1,nss
 126         call xdrfint(ixdrf, ihpb(j), iret)
 127         call xdrfint(ixdrf, jhpb(j), iret)
 128        enddo
 129        call xdrfint(ixdrf, nfrag, iret)
 130        call xdrfint(ixdrf, iset, iret)
 131        do i=1,nfrag
 132         call xdrffloat(ixdrf, qfrag(i), iret)
 133        enddo
 134        prec=10000.0
 135 c       print *,'nss=',nss,'nfrag=',nfrag,'iset=',iset
 136        isize=0
 137        call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
 138
 139
 140 c       write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 141 c       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
 142 c       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
 143 c       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
 144        if (mod(kk/ntraj,ifreq).eq.0 .and. mod(kk,ntraj).eq.itraj) then
 145          if (isize .ne. nres+nct-nnt+1) then
 146            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
 147          endif
 148          do i=1,nres
 149            do j=1,3
 150              c(j,i)=coord(j,i)
 151            enddo
 152          enddo
 153          ii = 0
 154          do i=nnt,nct
 155            ii = ii + 1
 156            do j=1,3
 157              c(j,i+nres)=coord(j,ii+nres)
 158            enddo
 159          enddo
 160 c Bring all chains to same box and to a common origin
 161          do ichain=1,nchain
 162 c           print *,"ichain",ichain
 163            do j=1,3
 164              cm(j,ichain)=0.0d0
 165            enddo
 166            do i=chain_border(1,ichain),chain_border(2,ichain)
 167 c             print *,i,c(:,i)
 168              do j=1,3
 169                cm(j,ichain)=cm(j,ichain)+c(j,i)
 170              enddo
 171            enddo
 172            do j=1,3
 173              cm(j,ichain)=cm(j,ichain)/chain_length(ichain)
 174            enddo
 175          enddo
 176 c         write (*,*) "CM"
 177 c         do i=1,nchain
 178 c           write (*,'(i5,3f10.5)') i,(cm(j,i),j=1,3)
 179 c         enddo
 180          do j=1,3
 181            lmin=0
 182            sumodl_min=1.0d10
 183            do i=0,3**nchain-1
 184              ii=i
 185              do ichain=1,nchain
 186                lshift(ichain)=mod(ii,3)-1
 187                ii=ii/3
 188              enddo
 189 c             print '(10i3)',(lshift(ichain),ichain=1,nchain)
 190              sumodl=0.0d0
 191              do ichain=1,nchain
 192                do jchain=1,ichain-1
 193                sumodl=sumodl+dabs(cm(j,ichain)+lshift(ichain)*boxsize(j)
 194      &             -cm(j,jchain)-lshift(jchain)*boxsize(j))
 195                enddo
 196              enddo
 197 c             print *,"sumodl",sumodl," sumodl_min",sumodl_min
 198              if (sumodl.lt.sumodl_min) then
 199                sumodl_min=sumodl
 200                lmin=lshift
 201              endif
 202            enddo
 203            do ichain=1,nchain
 204              cm(j,ichain)=cm(j,ichain)+lmin(ichain)*boxsize(j)
 205              shift_coord(j,ichain)=lmin(ichain)*boxsize(j)
 206            enddo
 207          enddo
 208          do ichain=1,nchain
 209            do j=1,3
 210              ccm(j)=ccm(j)+cm(j,ichain)*chain_length(ichain)
 211            enddo
 212          enddo
 213          do j=1,3
 214            ccm(j)=ccm(j)/nrestot
 215          enddo
 216 c         write (*,'(a,3f10.5)') "ccm     ",(ccm(j),j=1,3)
 217 c         write (*,'(a,3f10.5)') "ccm_prev",(ccm_prev(j),j=1,3)
 218 c         write (*,'(a,3f10.5)') "diff    ",(ccm(j)-ccm_prev(j),j=1,3)
 219 c         write (*,*) "shift_coord"
 220 c         do i=1,nchain
 221 c           write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
 222 c         enddo
 223 #ifdef PREVIOUS
 224          if (previous) then
 225            do j=1,3
 226              change=.true.
 227              do while (change)
 228              if (ccm(j)-ccm_prev(j).gt.boxsize(j)/2) then
 229                do ichain=1,nchain
 230                  shift_coord(j,ichain)=shift_coord(j,ichain)-boxsize(j)
 231                enddo
 232                ccm(j)=ccm(j)-boxsize(j)
 233              else if (ccm(j)-ccm_prev(j).lt.-boxsize(j)/2) then
 234                do ichain=1,nchain
 235                  shift_coord(j,ichain)=shift_coord(j,ichain)+boxsize(j)
 236                enddo
 237                ccm(j)=ccm(j)+boxsize(j)
 238              else
 239                change=.false.
 240              endif
 241            enddo
 242            enddo
 243            do j=1,3
 244              ccm_prev(j)=ccm(j)
 245            enddo
 246          else
 247            previous=.true.
 248            do j=1,3
 249              ccm_prev(j)=ccm(j)
 250            enddo
 251          endif
 252 #else
 253          do ichain=1,nchain
 254            do j=1,3
 255              shift_coord(j,ichain)=shift_coord(j,ichain)-ccm(j)
 256            enddo
 257          enddo
 258 #endif
 259 c         write (*,*) "shift_coord"
 260 c         do i=1,nchain
 261 c           write (*,'(i5,3f10.5)') i,(shift_coord(j,i),j=1,3)
 262 c         enddo
 263          do ichain=1,nchain
 264            do i=chain_border(1,ichain),chain_border(2,ichain)
 265              do j=1,3
 266                c(j,i)=c(j,i)+shift_coord(j,ichain)
 267                c(j,i+nres)=c(j,i+nres)+shift_coord(j,ichain)
 268              enddo
 269            enddo
 270          enddo
 271          etot=potE
 272          write (tytul,'(a1,i6,a8,f6.1,a6,f10.1)') "#",kk," t_bath ",
 273      &             t_bath, " time ",time
 274          call pdbout(etot,tytul,9)
 275        endif
 276       enddo
 277
 278       end