Fixed the "***" issue in PDB output format by shifting the center of each molecule...
[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       real*4 coord(3,2*maxres)
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
13       character*8 ucase
14       external ucase
15       logical oneletter
16       integer rescode
17       external rescode
18       logical iblnk
19       external iblnk
20       double precision cm(3)
21       
22       ifreq=1
23       is=1
24       ie=1000000000
25       if (iargc().lt.3) then
26         print '(2a)',
27      &    "Usage: xdrf2pdb one/three seqfile cxfile [freq]",
28      &    " [start] [end] [pdbfile]"
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,cfreq)
79         read (cfreq,*) is
80       endif
81       if (iargc().gt.5) then
82         call getarg(6,cfreq)
83         read (cfreq,*) ie
84       endif
85       if (iargc().gt.6) then
86         call getarg(7,pdbfile)
87       else
88         pdbfile=arg(:iext)//'.pdb'
89       endif
90       open(9,file=pdbfile)
91       nnt = 1
92       if (itype(1).eq.21) nnt = 2
93       nct=nres
94       if (itype(nres).eq.21) nct = nres-1
95       print *,"nnt",nnt," nct",nct
96       print *,"file",arg
97       call xdrfopen(ixdrf,arg, "r", iret)
98       if (iret.eq.0) stop
99       print *,"iret",iret
100       print *,"is",is," ie",ie
101       kk = 0
102       do while(is.eq.0 .or. kk.lt.ie) 
103        call xdrffloat(ixdrf, time, iret)
104        print *,"time",time," iret",iret
105        if(iret.eq.0) exit
106        kk = kk + 1
107        call xdrffloat(ixdrf, potE, iret)
108        call xdrffloat(ixdrf, uconst, iret)
109        print *,"potE",potE," uconst",uconst
110 #ifdef NEWUNRES
111        call xdrffloat(ixdrf, uconst_back, iret)
112        print *,"uconst_back",uconst_back
113 #endif
114        call xdrffloat(ixdrf, t_bath, iret)
115        print *,"t_bath",t_bath
116        call xdrfint(ixdrf, nss, iret) 
117        do j=1,nss
118         call xdrfint(ixdrf, ihpb(j), iret)
119         call xdrfint(ixdrf, jhpb(j), iret)
120        enddo
121        print *,"nss",nss
122        print *,"   ",(ihpb(j),"-",jhpb(j),j=1,nss)
123        call xdrfint(ixdrf, nfrag, iret)
124        do i=1,nfrag
125         call xdrffloat(ixdrf, qfrag(i), iret)
126        enddo
127        print *,"nfrag",nfrag
128        prec=10000.0
129
130        isize=0
131        call xdrf3dfcoord(ixdrf, coord, isize, prec, iret)
132
133
134 c       write (*,'(e15.8,2e15.5,f12.5,$)') time,potE,uconst,t_bath
135 c       write (*,'(i4,$)') nss,(ihpb(j),jhpb(j),j=1,nss)
136 c       write (*,'(i4,20f7.4)') nfrag,(qfrag(i),i=1,nfrag)
137 c       write (*,'(8f10.5)') ((coord(k,j),k=1,3),j=1,isize)
138        if (kk.ge.is .and. mod(kk,ifreq).eq.0) then
139          if (isize .ne. nres+nct-nnt+1) then
140            print *,"Error: inconsistent sizes",isize,nres+nct-nnt+1
141          endif
142          do i=1,nres
143            do j=1,3
144              c(j,i)=coord(j,i)
145            enddo
146          enddo
147          ii = 0
148          do i=nnt,nct
149            ii = ii + 1 
150            do j=1,3
151              c(j,i+nres)=coord(j,ii+nres)
152            enddo
153          enddo
154 c Calculate the CM
155          do j=1,3
156            cm(j)=0.0d0
157            do i=1,nres
158              cm(j)=cm(j)+c(j,i)
159            enddo
160            cm(j)=cm(j)/nres
161          enddo
162          do i=1,nres
163            do j=1,3
164              c(j,i)=c(j,i)-cm(j)
165              c(j,i+nres)=c(j,i+nres)-cm(j)
166            enddo
167          enddo
168          etot=potE
169          write (tytul,'(a,i6)') "Structure",kk
170          call pdbout(etot,tytul,9)
171        endif
172       enddo
173      
174       end