source/xdrfpdb/xdrf2pdb-m.F90

   1       program xdrfpdb_m
   2 !      implicit real*8 (a-h,o-z)
   3 !      include 'DIMENSIONS'
   4 !      include 'COMMON.CHAIN'
   5 !      include 'COMMON.INTERACT'
   6 !      include 'COMMON.SBRIDGE'
   7       use geometry_data, only: nres,c,nfrag,boxxsize,boxysize,boxzsize
   8       use energy_data, only: itype,nnt,nct,nss,ihpb,jhpb,uconst_back,molnum
   9       use io_base!, only: maxres,ucase,iblnk,rescode,pdbout
  10 !      use geometry, only: returnbox
  11       implicit none
  12       real(kind=4),allocatable,dimension(:,:) :: coord !(3,2*maxres)
  13       real(kind=4) :: prec,time,potE,uconst,t_bath,qfrag(100)
  14       real(kind=8) :: etot
  15       character(len=80) :: arg,seqfile,pdbfileX,boxtraj
  16       character(len=3),allocatable,dimension(:) :: sequenc !(maxres)
  17       character(len=50) :: tytul
  18       character(len=8) :: onethree,cfreq,cntraj,citraj
  19       character(len=3) :: licz
  20 !      external ucase
  21       logical oneletter !,iblnk
  22 !      integer rescode
  23 !      external rescode
  24       integer :: i,ii,j,k,kk,is,ie,ifreq,mnum,molec,iext,isize
  25       integer :: iret,ixdrf,itraj,ntraj
  26       print *,"begin program"
  27 !      if(.not.allocated(molnum)) allocate(molnum(maxres))
  28       ifreq=1
  29       nres=0
  30       mnum=1
  31       molec=1
  32       allocate(sequenc(maxres))
  33       allocate(molnum(maxres))
  34       allocate(itype(maxres,5))
  35       if (iargc().lt.6) then
  36         print '(2a)',&
  37          "Usage: xdrf2pdb-m one/three seqfile cxfile ntraj itraj",&
  38          "boxsize [pdbfile] [freq]"
  39         stop
  40       endif
  41       call getarg(1,onethree)
  42       onethree = ucase(onethree)
  43       if (onethree.eq.'ONE') then
  44         oneletter = .true.
  45       else if (onethree.eq.'THREE') then
  46         oneletter = .false.
  47       else
  48         print *,"ONE or THREE must be specified"
  49       endif
  50       call getarg(2,seqfile)
  51       open (1,file=seqfile,status='old')
  52       if (oneletter) then
  53         read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
  54    10   continue
  55         nres=i
  56         i=0
  57         do while (.not.iblnk(sequenc(i+1)(1:1)))
  58           i=i+1
  59         enddo
  60         nres=i
  61         do i=1,nres
  62           itype(i,mnum)=rescode(i,sequenc(i),1,molec)
  63         enddo
  64       else
  65         read(1,'(20(a3,i1))',end=11,err=11) (sequenc(i),molnum(i),i=1,maxres)
  66    11   continue
  67         nres=i
  68         i=0
  69         do while (.not.iblnk(sequenc(i+1)(1:1)))
  70           i=i+1
  71         enddo
  72         nres=i
  73         do i=1,nres
  74           itype(i,molnum(i))=rescode(i,sequenc(i),0,molnum(i))
  75         enddo
  76         print *,nres
  77         print '(20(a3,1x))',(sequenc(i),i=1,nres)
  78       endif
  79       call getarg(3,arg)
  80       iext = index(arg,'.cx') - 1
  81       if (iext.lt.0) then
  82         print *,"Error - not a cx file"
  83         stop
  84       endif
  85       call getarg(4,cntraj)
  86       read (cntraj,*) ntraj
  87       call getarg(5,citraj)
  88       read (citraj,*) itraj
  89       call getarg(6,boxtraj)
  90       read (boxtraj,*) boxxsize
  91       boxysize=boxxsize
  92       boxzsize=boxxsize
  93       if (iargc().gt.6) then
  94         call getarg(7,pdbfileX)
  95       else
  96         write(licz,'(bz,i3.3)') itraj
  97         pdbfileX=arg(:iext)//'_'//licz//'.pdb'
  98       endif
  99       if (iargc().gt.7) then
 100         call getarg(8,cfreq)
 101         read (cfreq,*) ifreq
 102       endif
 103 !      print *,"ifreq",ifreq," ntraj",ntraj," itraj",itraj
 104       open(9,file=pdbfileX)
 105       nnt = 1
 106       mnum=molnum(1)
 107       if (itype(1,mnum).eq.ntyp1) nnt = 2
 108       nct=nres
 109       mnum=molnum(nct)
 110       if (itype(nres,mnum).eq.ntyp1) nct = nres-1
 111       print *,"nnt",nnt," nct",nct
 112       call xdrfopen(ixdrf,arg, "r", iret)
 113       kk = 0
 114       allocate(coord(3,2*nres))
 115       allocate(c(3,2*nres))
 116       do while(.true.)
 117        call xdrffloat(ixdrf, time, iret)
 118        if(iret.eq.0) exit
 119        kk = kk + 1
 120        call xdrffloat(ixdrf, potE, iret)
 121        call xdrffloat(ixdrf, uconst, iret)
 122 #ifdef NEWUNRES
 123        call xdrffloat(ixdrf, uconst_back, iret)
 124 #endif
 125        call xdrffloat(ixdrf, t_bath, iret)
 126        call xdrfint(ixdrf, nss, iret)
 127        do j=1,nss
 128         call xdrfint(ixdrf, ihpb(j), iret)
 129         call xdrfint(ixdrf, jhpb(j), iret)
 130        enddo
 131        call xdrfint(ixdrf, nfrag, iret)
 132        do i=1,nfrag
 133         call xdrffloat(ixdrf, qfrag(i), iret)
 134        enddo
 135        prec=10000.0
 136
 137        isize=0
 138        call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
 139
 140
 141 !       write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 142 !       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
 143 !       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
 144        do j=1,isize
 145        write (*,'(3f10.5,i3)') (coord(k,j),k=1,3),j
 146        enddo
 147        if (mod(kk/ntraj,ifreq).eq.0 .and. mod(kk,ntraj).eq.itraj) then
 148          if (isize .ne. nres+nct-nnt+1) then
 149            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
 150          endif
 151          do i=1,nres
 152            do j=1,3
 153              c(j,i)=coord(j,i)
 154            enddo
 155          enddo
 156          ii = 0
 157          do i=nnt,nct
 158            ii = ii + 1
 159            do j=1,3
 160              c(j,i+nres)=coord(j,ii+nres)
 161            enddo
 162          enddo
 163          call returnbox
 164          etot=potE
 165          write (tytul,'(a,i6)') "Structure",kk
 166          call pdbout(etot,tytul,9)
 167        endif
 168       enddo
 169
 170       end
 171       subroutine returnbox
 172       use geometry_data, only: nres,c,nfrag,boxxsize,boxysize,boxzsize
 173       use energy_data, only: itype,nnt,nct,nss,ihpb,jhpb,uconst_back,molnum
 174       use io_base!, only: maxres,ucase,iblnk,rescode,pdbout
 175       integer :: allareout,i,j,k,nojumpval,chain_beg,mnum
 176       integer :: chain_end,ireturnval
 177       real*8 :: difference,xi,boxsize,x,xtemp,box2shift
 178
 179        do i=1,nres
 180          if (molnum(i).eq.5) then
 181           c(1,i)=dmod(c(1,i),boxxsize)
 182 !          if (c(1,i).lt.0) c(1,i)=c(1,i)+boxxsize
 183           x=c(1,i)-c(1,2)
 184           print *,"return box,before wrap",c(1,i)
 185          boxsize=boxxsize
 186            xtemp=dmod(x,boxsize)
 187       if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
 188         box2shift=xtemp-boxsize
 189       else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
 190         box2shift=xtemp+boxsize
 191       else
 192         box2shift=xtemp
 193       endif
 194           xi=box2shift
 195 !          print *,c(1,2),xi,xtemp,box2shift
 196           c(1,i)=c(1,2)+xi
 197           print *,"after",c(1,2),xi,xtemp,box2shift,c(1,i)
 198
 199           c(2,i)=dmod(c(2,i),boxysize)
 200 !          if (c(2,i).lt.0) c(2,i)=c(2,i)+boxysize
 201           x=c(2,i)-c(2,2)
 202          boxsize=boxysize
 203            xtemp=dmod(x,boxsize)
 204       if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
 205         box2shift=xtemp-boxsize
 206       else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
 207         box2shift=xtemp+boxsize
 208       else
 209         box2shift=xtemp
 210       endif
 211           xi=box2shift
 212           c(2,i)=c(2,2)+xi
 213 !         c(3,i)=dmod(c(3,i),boxzsize)
 214           if (c(3,i).lt.0) c(3,i)=c(3,i)+boxzsize
 215           x=c(3,i)-c(3,2)
 216          boxsize=boxzsize
 217            xtemp=dmod(x,boxsize)
 218       if (dabs(xtemp-boxsize).lt.dabs(xtemp)) then
 219         box2shift=xtemp-boxsize
 220       else if (dabs(xtemp+boxsize).lt.dabs(xtemp)) then
 221         box2shift=xtemp+boxsize
 222       else
 223         box2shift=xtemp
 224       endif
 225           xi=box2shift
 226           c(3,i)=c(3,2)+xi
 227
 228           c(1,i+nres)=dmod(c(1,i+nres),boxxsize)
 229           c(2,i+nres)=dmod(c(2,i+nres),boxysize)
 230           c(3,i+nres)=dmod(c(3,i+nres),boxzsize)
 231          endif
 232         enddo
 233         return
 234         end       subroutine returnbox
 235