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