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