source/xdrfpdb/xdrf2pdb.F90

   1       program xdrfpdb
   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
   8       use energy_data, only: itype,nnt,nct,nss,ihpb,jhpb,uconst_back
   9       use io_base!, only: maxres,ucase,iblnk,rescode,pdbout
  10
  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
  16       character(len=3),allocatable,dimension(:) :: sequenc !(maxres)
  17       character(len=50) tytul
  18       character(len=8) onethree,cfreq
  19 !      external ucase
  20       logical :: oneletter
  21 !      integer rescode
  22 !      external rescode
  23 !      logical iblnk
  24 !      external iblnk
  25 !el local variables
  26       integer :: i,ii,j,k,kk,is,ie,ifreq,mnum,molec,iext,isize
  27       integer :: iret,ixdrf
  28
  29       ifreq=1
  30       is=1
  31       ie=1000000000
  32       nres=0
  33       mnum=1
  34       molec=1
  35       allocate(sequenc(maxres))
  36       allocate(itype(maxres,5))
  37 !      allocate(c(3,2*maxres+2))
  38       if (iargc().lt.3) then
  39         print '(2a)',&
  40           "Usage: xdrf2pdb one/three seqfile cxfile [freq]",&
  41           " [start] [end] [pdbfile]"
  42         stop
  43       endif
  44       call getarg(1,onethree)
  45       onethree = ucase(onethree)
  46       if (onethree.eq.'ONE') then
  47         oneletter = .true.
  48       else if (onethree.eq.'THREE') then
  49         oneletter = .false.
  50       else
  51         print *,"ONE or THREE must be specified"
  52       endif
  53       call getarg(2,seqfile)
  54       open (1,file=seqfile,status='old')
  55       if (oneletter) then
  56         read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
  57    10   continue
  58         nres=i
  59         i=0
  60         do while (.not.iblnk(sequenc(i+1)(1:1)))
  61           i=i+1
  62         enddo
  63         nres=i
  64         do i=1,nres
  65           itype(i,mnum)=rescode(i,sequenc(i),1,molec)
  66         enddo
  67       else
  68         read(1,'(20(a3,1x))',end=11,err=11) (sequenc(i),i=1,maxres)
  69    11   continue
  70         nres=i
  71         i=0
  72         do while (.not.iblnk(sequenc(i+1)(1:1)))
  73           i=i+1
  74         enddo
  75         nres=i
  76         do i=1,nres
  77           itype(i,mnum)=rescode(i,sequenc(i),0,molec)
  78         enddo
  79       endif
  80       call getarg(3,arg)
  81       iext = index(arg,'.cx') - 1
  82       if (iext.lt.0) then
  83         print *,"Error - not a cx file"
  84         stop
  85       endif
  86       if (iargc().gt.3) then
  87         call getarg(4,cfreq)
  88         read (cfreq,*) ifreq
  89       endif
  90       if (iargc().gt.4) then
  91         call getarg(5,cfreq)
  92         read (cfreq,*) is
  93       endif
  94       if (iargc().gt.5) then
  95         call getarg(6,cfreq)
  96         read (cfreq,*) ie
  97       endif
  98       if (iargc().gt.6) then
  99         call getarg(7,pdbfileX)
 100       else
 101         pdbfileX=arg(:iext)//'.pdb'
 102       endif
 103       open(9,file=pdbfileX)
 104       nnt = 1
 105       if (itype(1,mnum).eq.ntyp1) nnt = 2
 106       nct=nres
 107       if (itype(nres,mnum).eq.ntyp1) nct = nres-1
 108       print *,"nnt",nnt," nct",nct," nres",nres
 109       print *,"file",arg
 110       call xdrfopen(ixdrf,arg, "r", iret)
 111       if (iret.eq.0) stop
 112       print *,"iret",iret
 113       print *,"is",is," ie",ie
 114       kk = 0
 115       allocate(coord(3,2*nres))
 116       allocate(c(3,2*nres))
 117       do while(is.eq.0 .or. kk.lt.ie)
 118        call xdrffloat(ixdrf, time, iret)
 119        print *,"time",time," iret",iret
 120        if(iret.eq.0) exit
 121        kk = kk + 1
 122        call xdrffloat(ixdrf, potE, iret)
 123        call xdrffloat(ixdrf, uconst, iret)
 124        print *,"potE",potE," uconst",uconst
 125        call xdrffloat(ixdrf, uconst_back, iret)
 126        call xdrffloat(ixdrf, t_bath, iret)
 127        print *,"t_bath",t_bath
 128 !#ifdef NEWUNRES
 129 !#endif
 130        print *,"uconst_back",uconst_back
 131 !       call xdrffloat(ixdrf, t_bath, iret)
 132 !       print *,"t_bath",t_bath
 133        call xdrfint(ixdrf, nss, iret)
 134        do j=1,nss
 135         call xdrfint(ixdrf, ihpb(j), iret)
 136         call xdrfint(ixdrf, jhpb(j), iret)
 137        enddo
 138        print *,"nss",nss
 139 !       call xdrffloat(ixdrf, t_bath, iret)
 140 !       print *,"t_bath",t_bath
 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
 151
 152 !       write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 153 !       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
 154 !       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
 155 !       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
 156        if (kk.ge.is .and. mod(kk,ifreq).eq.0) then
 157          if (isize .ne. nres+nct-nnt+1) then
 158            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
 159          endif
 160          do i=1,nres
 161            do j=1,3
 162              c(j,i)=coord(j,i)
 163            enddo
 164          enddo
 165          ii = 0
 166          do i=nnt,nct
 167            ii = ii + 1
 168            do j=1,3
 169              c(j,i+nres)=coord(j,ii+nres)
 170            enddo
 171          enddo
 172 !       do j=1,nres
 173 !         write (*,'(8f10.5)') (c(k,j),k=1,3),(c(i,j+nres),i=1,3)
 174 !       enddo
 175          etot=potE
 176          write (tytul,'(a,i6,a,f8.3)') "Structure",kk,"Temp",t_bath
 177          call pdbout(etot,tytul,9)
 178        endif
 179       enddo
 180
 181       end