Merge branch 'devel' of mmka.chem.univ.gda.pl:unres into devel
[unres.git] / source / xdrfpdb / src / xdrf2ang.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,2*maxres)
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 <<<<<<< HEAD
18       logical oneletter,iblnk
19 =======
20       logical oneletter, iblnk
21 >>>>>>> aebadf9023e437f497f92dfcf2303c16f60917c1
22       integer rescode
23       external rescode
24       
25       pi = 4.0d0*datan(1.0d0)
26       deg2rad=pi/180.0d0
27       rad2deg=1.0d0/deg2rad
28
29       ifreq=1
30       is=1
31       ie=1000000000
32       if (iargc().lt.3) then
33         print '(2a)',
34      &    "Usage: xdrf2ang one/three seqfile cxfile [freq]",
35      &    " [start] [end] [angfile]"
36         stop
37       endif
38       call getarg(1,onethree)
39       onethree = ucase(onethree)
40       if (onethree.eq.'ONE') then
41         oneletter = .true.
42       else if (onethree.eq.'THREE') then
43         oneletter = .false.
44       else
45         print *,"ONE or THREE must be specified"
46       endif
47       call getarg(2,seqfile)
48       open (1,file=seqfile,status='old')
49       if (oneletter) then
50         read(1,'(80a1)',end=10,err=10) (sequenc(i)(1:1),i=1,maxres)
51    10   continue
52         nres=i
53         i=0
54         do while (.not.iblnk(sequenc(i+1)(1:1)))
55 c        do while (.not.(iblnk(sequenc(i+1)(1:1)) == 0) )
56           i=i+1
57         enddo 
58         nres=i
59         do i=1,nres
60           itype(i)=rescode(i,sequenc(i),1)
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 c        do while (.not.(iblnk(sequenc(i+1)(1:1)) == 0) )
69           i=i+1
70         enddo 
71         nres=i
72         do i=1,nres
73           itype(i)=rescode(i,sequenc(i),0)
74         enddo
75       endif
76       call getarg(3,arg)
77       iext = index(arg,'.cx') - 1
78       if (iext.lt.0) then
79         print *,"Error - not a cx file"
80         stop
81       endif
82       if (iargc().gt.3) then
83         call getarg(4,cfreq)
84         read (cfreq,*) ifreq
85       endif
86       if (iargc().gt.4) then
87         call getarg(5,cfreq)
88         read (cfreq,*) is
89       endif
90       if (iargc().gt.5) then
91         call getarg(6,cfreq)
92         read (cfreq,*) ie
93       endif
94       if (iargc().gt.6) then
95         call getarg(7,angfile)
96       else
97         angfile=arg(:iext)//'.ang'
98       endif
99       open(9,file=angfile)
100       nnt = 1
101       if (itype(1).eq.21) nnt = 2
102       nct=nres
103       if (itype(nres).eq.21) nct = nres-1
104       print *,"nnt",nnt," nct",nct
105       print *,"file",arg
106       call xdrfopen(ixdrf,arg, "r", iret)
107       if (iret.eq.0) stop
108       print *,"iret",iret
109       print *,"is",is," ie",ie
110       kk = 0
111       do while(is.eq.0 .or. kk.lt.ie) 
112 c       call xdrfint(ixdrf, ic, iret)
113 c       print *,ic
114        call xdrffloat(ixdrf, time, iret)
115        print *,"time",time," iret",iret
116        if(iret.eq.0) exit
117        kk = kk + 1
118        call xdrffloat(ixdrf, potE, iret)
119        call xdrffloat(ixdrf, uconst, iret)
120        call xdrffloat(ixdrf, uconst_back, iret)
121 c       print *,"potE",potE," uconst",uconst," uconst_back",uconst_back
122        call xdrffloat(ixdrf, t_bath, iret)
123 c       print *,"t_bath",t_bath
124        call xdrfint(ixdrf, nss, iret) 
125 c       print *,"nss",nss
126        do j=1,nss
127         call xdrfint(ixdrf, ihpb(j), iret)
128         call xdrfint(ixdrf, jhpb(j), iret)
129        enddo
130        call xdrfint(ixdrf, nfrag, iret)
131 c       print *,"nfrag",nfrag
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 c       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
142 c       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
143 c       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
144        if (mod(kk,ifreq).eq.0) then
145          if (isize .ne. nres+nct-nnt+1) then
146            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
147          endif
148          write (9,'(i10,e15.8,2e15.5,f12.5)') kk,time,potE,uconst,t_bath
149          do i=1,nres
150            do j=1,3
151              c(j,i)=coord(j,i)
152            enddo
153          enddo
154          ii = 0
155          do i=nnt,nct
156            ii = ii + 1 
157            do j=1,3
158              c(j,i+nres)=coord(j,ii+nres)
159            enddo
160          enddo
161          etot=potE
162          do i=2,nres-2
163            write (9,'(a3,1x,2f10.3)')
164      &       restyp(itype(i)),alpha(i-1,i,i+1)*rad2deg,
165      &       beta(i-1,i,i+1,i+2)*rad2deg
166          enddo
167          write (9,'(a3,1x,f10.3)')
168      &       restyp(itype(nres-1)),alpha(nres-2,nres-1,nres)*rad2deg
169        endif
170       enddo
171      
172       end