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
   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,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
  27       ifreq=1
  28       nres=0
  29       mnum=1
  30       molec=1
  31       allocate(sequenc(maxres))
  32       allocate(itype(maxres,5))
  33       if (iargc().lt.5) then
  34         print '(2a)',&
  35          "Usage: xdrf2pdb-m one/three seqfile cxfile ntraj itraj",&
  36          " [pdbfile] [freq]"
  37         stop
  38       endif
  39       call getarg(1,onethree)
  40       onethree = ucase(onethree)
  41       if (onethree.eq.'ONE') then
  42         oneletter = .true.
  43       else if (onethree.eq.'THREE') then
  44         oneletter = .false.
  45       else
  46         print *,"ONE or THREE must be specified"
  47       endif
  48       call getarg(2,seqfile)
  49       open (1,file=seqfile,status='old')
  50       if (oneletter) then
  51         read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
  52    10   continue
  53         nres=i
  54         i=0
  55         do while (.not.iblnk(sequenc(i+1)(1:1)))
  56           i=i+1
  57         enddo
  58         nres=i
  59         do i=1,nres
  60           itype(i,mnum)=rescode(i,sequenc(i),1,molec)
  61         enddo
  62       else
  63         read(1,'(20(a3,1x))',end=11,err=11) (sequenc(i),i=1,maxres)
  64    11   continue
  65         nres=i
  66         i=0
  67         do while (.not.iblnk(sequenc(i+1)(1:1)))
  68           i=i+1
  69         enddo
  70         nres=i
  71         do i=1,nres
  72           itype(i,mnum)=rescode(i,sequenc(i),0,molec)
  73         enddo
  74         print *,nres
  75         print '(20(a3,1x))',(sequenc(i),i=1,nres)
  76       endif
  77       call getarg(3,arg)
  78       iext = index(arg,'.cx') - 1
  79       if (iext.lt.0) then
  80         print *,"Error - not a cx file"
  81         stop
  82       endif
  83       call getarg(4,cntraj)
  84       read (cntraj,*) ntraj
  85       call getarg(5,citraj)
  86       read (citraj,*) itraj
  87       if (iargc().gt.5) then
  88         call getarg(6,pdbfileX)
  89       else
  90         write(licz,'(bz,i3.3)') itraj
  91         pdbfileX=arg(:iext)//'_'//licz//'.pdb'
  92       endif
  93       if (iargc().gt.6) then
  94         call getarg(7,cfreq)
  95         read (cfreq,*) ifreq
  96       endif
  97 !      print *,"ifreq",ifreq," ntraj",ntraj," itraj",itraj
  98       open(9,file=pdbfileX)
  99       nnt = 1
 100       if (itype(1,mnum).eq.ntyp1) nnt = 2
 101       nct=nres
 102       if (itype(nres,mnum).eq.ntyp1) nct = nres-1
 103       print *,"nnt",nnt," nct",nct
 104       call xdrfopen(ixdrf,arg, "r", iret)
 105       kk = 0
 106       allocate(coord(3,2*nres))
 107       allocate(c(3,2*nres))
 108       do while(.true.)
 109        call xdrffloat(ixdrf, time, iret)
 110        if(iret.eq.0) exit
 111        kk = kk + 1
 112        call xdrffloat(ixdrf, potE, iret)
 113        call xdrffloat(ixdrf, uconst, iret)
 114 #ifdef NEWUNRES
 115        call xdrffloat(ixdrf, uconst_back, iret)
 116 #endif
 117        call xdrffloat(ixdrf, t_bath, iret)
 118        call xdrfint(ixdrf, nss, iret)
 119        do j=1,nss
 120         call xdrfint(ixdrf, ihpb(j), iret)
 121         call xdrfint(ixdrf, jhpb(j), iret)
 122        enddo
 123        call xdrfint(ixdrf, nfrag, iret)
 124        do i=1,nfrag
 125         call xdrffloat(ixdrf, qfrag(i), iret)
 126        enddo
 127        prec=10000.0
 128
 129        isize=0
 130        call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
 131
 132
 133 !       write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
 134 !       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
 135 !       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
 136 !       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
 137        if (mod(kk/ntraj,ifreq).eq.0 .and. mod(kk,ntraj).eq.itraj) then
 138          if (isize .ne. nres+nct-nnt+1) then
 139            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
 140          endif
 141          do i=1,nres
 142            do j=1,3
 143              c(j,i)=coord(j,i)
 144            enddo
 145          enddo
 146          ii = 0
 147          do i=nnt,nct
 148            ii = ii + 1
 149            do j=1,3
 150              c(j,i+nres)=coord(j,ii+nres)
 151            enddo
 152          enddo
 153          etot=potE
 154          write (tytul,'(a,i6)') "Structure",kk
 155          call pdbout(etot,tytul,9)
 156        endif
 157       enddo
 158
 159       end