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